src/SMDS/SMDS_VolumeTool.cxx

   1 // Copyright (C) 2007-2023  CEA, EDF, 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_Mesh.hxx"
  36
  37 #include <utilities.h>
  38
  39 #include <map>
  40 #include <limits>
  41 #include <cmath>
  42 #include <cstring>
  43 #include <numeric>
  44 #include <algorithm>
  45 #include <array>
  46
  47 namespace
  48 {
  49 // ======================================================
  50 // Node indices in faces depending on volume orientation
  51 // making most faces normals external
  52 // ======================================================
  53 // For all elements, 0-th face is bottom based on the first nodes.
  54 // For prismatic elements (tetra,hexa,prisms), 1-th face is a top one.
  55 // For all elements, side faces follow order of bottom nodes
  56 // ======================================================
  57
  58 /*
  59 //           N3
  60 //           +
  61 //          /|\
  62 //         / | \
  63 //        /  |  \
  64 //    N0 +---|---+ N1                TETRAHEDRON
  65 //       \   |   /
  66 //        \  |  /
  67 //         \ | /
  68 //          \|/
  69 //           +
  70 //           N2
  71 */
  72 static int Tetra_F [4][4] = { // FORWARD == EXTERNAL
  73   { 0, 1, 2, 0 },              // All faces have external normals
  74   { 0, 3, 1, 0 },
  75   { 1, 3, 2, 1 },
  76   { 0, 2, 3, 0 }};
  77 static int Tetra_RE [4][4] = { // REVERSED -> FORWARD (EXTERNAL)
  78   { 0, 2, 1, 0 },              // All faces have external normals
  79   { 0, 1, 3, 0 },
  80   { 1, 2, 3, 1 },
  81   { 0, 3, 2, 0 }};
  82 static int Tetra_nbN [] = { 3, 3, 3, 3 };
  83
  84 /*
  85 //    N1 +---------+ N2
  86 //       | \     / |
  87 //       |  \   /  |
  88 //       |   \ /   |
  89 //       |   /+\   |     PYRAMID
  90 //       |  / N4\  |
  91 //       | /     \ |
  92 //       |/       \|
  93 //    N0 +---------+ N3
  94 */
  95 static int Pyramid_F [5][5] = { // FORWARD == EXTERNAL
  96   { 0, 1, 2, 3, 0 },            // All faces have external normals
  97   { 0, 4, 1, 0, 4 },
  98   { 1, 4, 2, 1, 4 },
  99   { 2, 4, 3, 2, 4 },
 100   { 3, 4, 0, 3, 4 }
 101 };
 102 static int Pyramid_RE [5][5] = { // REVERSED -> FORWARD (EXTERNAL)
 103   { 0, 3, 2, 1, 0 },             // All faces but a bottom have external normals
 104   { 0, 1, 4, 0, 4 },
 105   { 1, 2, 4, 1, 4 },
 106   { 2, 3, 4, 2, 4 },
 107   { 3, 0, 4, 3, 4 }};
 108 static int Pyramid_nbN [] = { 4, 3, 3, 3, 3 };
 109
 110 /*
 111 //            + N4
 112 //           /|\
 113 //          / | \
 114 //         /  |  \
 115 //        /   |   \
 116 //    N3 +---------+ N5
 117 //       |    |    |
 118 //       |    + N1 |
 119 //       |   / \   |                PENTAHEDRON
 120 //       |  /   \  |
 121 //       | /     \ |
 122 //       |/       \|
 123 //    N0 +---------+ N2
 124 */
 125 static int Penta_F [5][5] = { // FORWARD
 126   { 0, 1, 2, 0, 0 },          // All faces have external normals
 127   { 3, 5, 4, 3, 3 },          // 0 is bottom, 1 is top face
 128   { 0, 3, 4, 1, 0 },
 129   { 1, 4, 5, 2, 1 },
 130   { 0, 2, 5, 3, 0 }};
 131 static int Penta_RE [5][5] = { // REVERSED -> EXTERNAL
 132   { 0, 2, 1, 0, 0 },
 133   { 3, 4, 5, 3, 3 },
 134   { 0, 1, 4, 3, 0 },
 135   { 1, 2, 5, 4, 1 },
 136   { 0, 3, 5, 2, 0 }};
 137 static int Penta_nbN [] = { 3, 3, 4, 4, 4 };
 138
 139 /*
 140 //         N5+----------+N6
 141 //          /|         /|
 142 //         / |        / |
 143 //        /  |       /  |
 144 //     N4+----------+N7 |
 145 //       |   |      |   |           HEXAHEDRON
 146 //       | N1+------|---+N2
 147 //       |  /       |  /
 148 //       | /        | /
 149 //       |/         |/
 150 //     N0+----------+N3
 151 */
 152 static int Hexa_F [6][5] = { // FORWARD
 153   { 0, 1, 2, 3, 0 },
 154   { 4, 7, 6, 5, 4 },          // all face normals are external
 155   { 0, 4, 5, 1, 0 },
 156   { 1, 5, 6, 2, 1 },
 157   { 3, 2, 6, 7, 3 },
 158   { 0, 3, 7, 4, 0 }};
 159 static int Hexa_RE [6][5] = { // REVERSED -> EXTERNAL
 160   { 0, 3, 2, 1, 0 },
 161   { 4, 5, 6, 7, 4 },          // all face normals are external
 162   { 0, 1, 5, 4, 0 },
 163   { 1, 2, 6, 5, 1 },
 164   { 3, 7, 6, 2, 3 },
 165   { 0, 4, 7, 3, 0 }};
 166 static int Hexa_nbN [] = { 4, 4, 4, 4, 4, 4 };
 167 static int Hexa_oppF[] = { 1, 0, 4, 5, 2, 3 }; // oppopsite facet indices
 168
 169 /*
 170 //      N8 +------+ N9
 171 //        /        \
 172 //       /          \
 173 //   N7 +            + N10
 174 //       \          /
 175 //        \        /
 176 //      N6 +------+ N11
 177 //                             HEXAGONAL PRISM
 178 //      N2 +------+ N3
 179 //        /        \
 180 //       /          \
 181 //   N1 +            + N4
 182 //       \          /
 183 //        \        /
 184 //      N0 +------+ N5
 185 */
 186 static int HexPrism_F [8][7] = { // FORWARD
 187   { 0, 1, 2, 3, 4, 5, 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_RE [8][7] = { // REVERSED -> EXTERNAL
 196   { 0, 5, 4, 3, 2, 1, 0 },
 197   { 6,11,10, 9, 8, 7, 6 },
 198   { 0, 6, 7, 1, 0, 0, 0 },
 199   { 1, 7, 8, 2, 1, 1, 1 },
 200   { 2, 8, 9, 3, 2, 2, 2 },
 201   { 3, 9,10, 4, 3, 3, 3 },
 202   { 4,10,11, 5, 4, 4, 4 },
 203   { 5,11, 6, 0, 5, 5, 5 }};
 204 static int HexPrism_nbN [] = { 6, 6, 4, 4, 4, 4, 4, 4 };
 205
 206
 207 /*
 208 //           N3
 209 //           +
 210 //          /|\
 211 //        7/ | \8
 212 //        /  |4 \                    QUADRATIC
 213 //    N0 +---|---+ N1                TETRAHEDRON
 214 //       \   +9  /
 215 //        \  |  /
 216 //        6\ | /5
 217 //          \|/
 218 //           +
 219 //           N2
 220 */
 221 static int QuadTetra_F [4][7] = { // FORWARD
 222   { 0, 4, 1, 5, 2, 6, 0 },        // All faces have external normals
 223   { 0, 7, 3, 8, 1, 4, 0 },
 224   { 1, 8, 3, 9, 2, 5, 1 },
 225   { 0, 6, 2, 9, 3, 7, 0 }};
 226 static int QuadTetra_RE [4][7] = { // REVERSED -> FORWARD (EXTERNAL)
 227   { 0, 6, 2, 5, 1, 4, 0 },         // All faces have external normals
 228   { 0, 4, 1, 8, 3, 7, 0 },
 229   { 1, 5, 2, 9, 3, 8, 1 },
 230   { 0, 7, 3, 9, 2, 6, 0 }};
 231 static int QuadTetra_nbN [] = { 6, 6, 6, 6 };
 232
 233 //
 234 //     QUADRATIC
 235 //     PYRAMID
 236 //
 237 //            +4
 238 //
 239 //
 240 //       10+-----+11
 241 //         |     |        9 - middle point for (0,4) etc.
 242 //         |     |
 243 //        9+-----+12
 244 //
 245 //            6
 246 //      1+----+----+2
 247 //       |         |
 248 //       |         |
 249 //      5+         +7
 250 //       |         |
 251 //       |         |
 252 //      0+----+----+3
 253 //            8
 254 static int QuadPyram_F [5][9] = {  // FORWARD
 255   { 0, 5, 1, 6, 2, 7, 3, 8, 0 },   // All faces have external normals
 256   { 0, 9, 4, 10,1, 5, 0, 4, 4 },
 257   { 1, 10,4, 11,2, 6, 1, 4, 4 },
 258   { 2, 11,4, 12,3, 7, 2, 4, 4 },
 259   { 3, 12,4, 9, 0, 8, 3, 4, 4 }};
 260 static int QuadPyram_RE [5][9] = { // REVERSED -> FORWARD (EXTERNAL)
 261   { 0, 8, 3, 7, 2, 6, 1, 5, 0 },   // All faces but a bottom have external normals
 262   { 0, 5, 1, 10,4, 9, 0, 4, 4 },
 263   { 1, 6, 2, 11,4, 10,1, 4, 4 },
 264   { 2, 7, 3, 12,4, 11,2, 4, 4 },
 265   { 3, 8, 0, 9, 4, 12,3, 4, 4 }};
 266 static int QuadPyram_nbN [] = { 8, 6, 6, 6, 6 };
 267
 268 /*
 269 //            + N4
 270 //           /|\
 271 //         9/ | \10
 272 //         /  |  \
 273 //        /   |   \
 274 //    N3 +----+----+ N5
 275 //       |    |11  |
 276 //       |    |    |
 277 //       |    +13  |                QUADRATIC
 278 //       |    |    |                PENTAHEDRON
 279 //     12+    |    +14
 280 //       |    |    |
 281 //       |    |    |
 282 //       |    + N1 |
 283 //       |   / \   |
 284 //       | 6/   \7 |
 285 //       | /     \ |
 286 //       |/       \|
 287 //    N0 +---------+ N2
 288 //            8
 289 */
 290 static int QuadPenta_F [5][9] = {  // FORWARD
 291   { 0, 6, 1, 7, 2, 8, 0, 0, 0 },
 292   { 3, 11,5, 10,4, 9, 3, 3, 3 },
 293   { 0, 12,3, 9, 4, 13,1, 6, 0 },
 294   { 1, 13,4, 10,5, 14,2, 7, 1 },
 295   { 0, 8, 2, 14,5, 11,3, 12,0 }};
 296 static int QuadPenta_RE [5][9] = { // REVERSED -> EXTERNAL
 297   { 0, 8, 2, 7, 1, 6, 0, 0, 0 },
 298   { 3, 9, 4, 10,5, 11,3, 3, 3 },
 299   { 0, 6, 1, 13,4, 9, 3, 12,0 },
 300   { 1, 7, 2, 14,5, 10,4, 13,1 },
 301   { 0, 12,3, 11,5, 14,2, 8, 0 }};
 302 static int QuadPenta_nbN [] = { 6, 6, 8, 8, 8 };
 303
 304 /*
 305 //                 13
 306 //         N5+-----+-----+N6                          +-----+-----+
 307 //          /|          /|                           /|          /|
 308 //       12+ |       14+ |                          + |   +25   + |
 309 //        /  |        /  |                         /  |        /  |
 310 //     N4+-----+-----+N7 |       QUADRATIC        +-----+-----+   |  Central nodes
 311 //       |   | 15    |   |       HEXAHEDRON       |   |       |   |  of tri-quadratic
 312 //       |   |       |   |                        |   |       |   |  HEXAHEDRON
 313 //       | 17+       |   +18                      |   +   22+ |   +
 314 //       |   |       |   |                        |21 |       |   |
 315 //       |   |       |   |                        | + | 26+   | + |
 316 //       |   |       |   |                        |   |       |23 |
 317 //     16+   |       +19 |                        +   | +24   +   |
 318 //       |   |       |   |                        |   |       |   |
 319 //       |   |     9 |   |                        |   |       |   |
 320 //       | N1+-----+-|---+N2                      |   +-----+-|---+
 321 //       |  /        |  /                         |  /        |  /
 322 //       | +8        | +10                        | +   20+   | +
 323 //       |/          |/                           |/          |/
 324 //     N0+-----+-----+N3                          +-----+-----+
 325 //             11
 326 */
 327 static int QuadHexa_F [6][9] = {  // FORWARD
 328   { 0, 8, 1, 9, 2, 10,3, 11,0 },   // all face normals are external,
 329   { 4, 15,7, 14,6, 13,5, 12,4 },
 330   { 0, 16,4, 12,5, 17,1, 8, 0 },
 331   { 1, 17,5, 13,6, 18,2, 9, 1 },
 332   { 3, 10,2, 18,6, 14,7, 19,3 },
 333   { 0, 11,3, 19,7, 15,4, 16,0 }};
 334 static int QuadHexa_RE [6][9] = {  // REVERSED -> EXTERNAL
 335   { 0, 11,3, 10,2, 9, 1, 8, 0 },   // all face normals are external
 336   { 4, 12,5, 13,6, 14,7, 15,4 },
 337   { 0, 8, 1, 17,5, 12,4, 16,0 },
 338   { 1, 9, 2, 18,6, 13,5, 17,1 },
 339   { 3, 19,7, 14,6, 18,2, 10,3 },
 340   { 0, 16,4, 15,7, 19,3, 11,0 }};
 341 static int QuadHexa_nbN [] = { 8, 8, 8, 8, 8, 8 };
 342
 343 static int TriQuadHexa_F [6][9] = {  // FORWARD
 344   { 0, 8, 1, 9, 2, 10,3, 11, 20 },   // all face normals are external
 345   { 4, 15,7, 14,6, 13,5, 12, 25 },
 346   { 0, 16,4, 12,5, 17,1, 8,  21 },
 347   { 1, 17,5, 13,6, 18,2, 9,  22 },
 348   { 3, 10,2, 18,6, 14,7, 19, 23 },
 349   { 0, 11,3, 19,7, 15,4, 16, 24 }};
 350 static int TriQuadHexa_RE [6][9] = {  // REVERSED -> EXTERNAL
 351   { 0, 11,3, 10,2, 9, 1, 8,  20 },   // opposite faces are neighbouring,
 352   { 4, 12,5, 13,6, 14,7, 15, 25 },   // all face normals are external
 353   { 0, 8, 1, 17,5, 12,4, 16, 21 },
 354   { 1, 9, 2, 18,6, 13,5, 17, 22 },
 355   { 3, 19,7, 14,6, 18,2, 10, 23 },
 356   { 0, 16,4, 15,7, 19,3, 11, 24 }};
 357 static int TriQuadHexa_nbN [] = { 9, 9, 9, 9, 9, 9 };
 358
 359
 360 // ========================================================
 361 // to perform some calculations without linkage to CASCADE
 362 // ========================================================
 363 struct XYZ {
 364   double x;
 365   double y;
 366   double z;
 367   XYZ()                               { x = 0; y = 0; z = 0; }
 368   XYZ( double X, double Y, double Z ) { x = X; y = Y; z = Z; }
 369   XYZ( const XYZ& other )             { x = other.x; y = other.y; z = other.z; }
 370   XYZ( const SMDS_MeshNode* n )       { x = n->X(); y = n->Y(); z = n->Z(); }
 371   double* data()                      { return &x; }
 372   inline XYZ operator-( const XYZ& other );
 373   inline XYZ operator+( const XYZ& other );
 374   inline XYZ Crossed( const XYZ& other );
 375   inline XYZ operator-();
 376   inline double Dot( const XYZ& other );
 377   inline double Magnitude();
 378   inline double SquareMagnitude();
 379   inline XYZ Normalize();
 380 };
 381 inline XYZ XYZ::operator-( const XYZ& Right ) {
 382   return XYZ(x - Right.x, y - Right.y, z - Right.z);
 383 }
 384 inline XYZ XYZ::operator-() {
 385   return XYZ(-x,-y,-z);
 386 }
 387 inline XYZ XYZ::operator+( const XYZ& Right ) {
 388   return XYZ(x + Right.x, y + Right.y, z + Right.z);
 389 }
 390 inline XYZ XYZ::Crossed( const XYZ& Right ) {
 391   return XYZ (y * Right.z - z * Right.y,
 392               z * Right.x - x * Right.z,
 393               x * Right.y - y * Right.x);
 394 }
 395 inline double XYZ::Dot( const XYZ& Other ) {
 396   return(x * Other.x + y * Other.y + z * Other.z);
 397 }
 398 inline double XYZ::Magnitude() {
 399   return sqrt (x * x + y * y + z * z);
 400 }
 401 inline double XYZ::SquareMagnitude() {
 402   return (x * x + y * y + z * z);
 403 }
 404 inline XYZ XYZ::Normalize() {
 405   double magnitude = Magnitude();
 406   if ( magnitude != 0.0 )
 407     return XYZ(x /= magnitude,y /= magnitude,z /= magnitude );
 408   else
 409     return XYZ(x,y,z);
 410 }
 411
 412   //================================================================================
 413   /*!
 414    * \brief Return linear type corresponding to a quadratic one
 415    */
 416   //================================================================================
 417
 418   SMDS_VolumeTool::VolumeType quadToLinear(SMDS_VolumeTool::VolumeType quadType)
 419   {
 420     SMDS_VolumeTool::VolumeType linType = SMDS_VolumeTool::VolumeType( int(quadType)-4 );
 421     const int           nbCornersByQuad = SMDS_VolumeTool::NbCornerNodes( quadType );
 422     if ( SMDS_VolumeTool::NbCornerNodes( linType ) == nbCornersByQuad )
 423       return linType;
 424
 425     int iLin = 0;
 426     for ( ; iLin < SMDS_VolumeTool::NB_VOLUME_TYPES; ++iLin )
 427       if ( SMDS_VolumeTool::NbCornerNodes( SMDS_VolumeTool::VolumeType( iLin )) == nbCornersByQuad)
 428         return SMDS_VolumeTool::VolumeType( iLin );
 429
 430     return SMDS_VolumeTool::UNKNOWN;
 431   }
 432
 433 } // namespace
 434
 435 //================================================================================
 436 /*!
 437  * \brief Saver/restorer of a SMDS_VolumeTool::myCurFace
 438  */
 439 //================================================================================
 440
 441 struct SMDS_VolumeTool::SaveFacet
 442 {
 443   SMDS_VolumeTool::Facet  mySaved;
 444   SMDS_VolumeTool::Facet& myToRestore;
 445   SaveFacet( SMDS_VolumeTool::Facet& facet ): myToRestore( facet )
 446   {
 447     mySaved = facet;
 448     mySaved.myNodes.swap( facet.myNodes );
 449   }
 450   ~SaveFacet()
 451   {
 452     if ( myToRestore.myIndex != mySaved.myIndex )
 453       myToRestore = mySaved;
 454     myToRestore.myNodes.swap( mySaved.myNodes );
 455   }
 456 };
 457
 458 //=======================================================================
 459 //function : SMDS_VolumeTool
 460 //purpose  :
 461 //=======================================================================
 462
 463 SMDS_VolumeTool::SMDS_VolumeTool ()
 464 {
 465   Set( 0 );
 466 }
 467
 468 //=======================================================================
 469 //function : SMDS_VolumeTool
 470 //purpose  :
 471 //=======================================================================
 472
 473 SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume,
 474                                   const bool              ignoreCentralNodes)
 475 {
 476   Set( theVolume, ignoreCentralNodes );
 477 }
 478
 479 //=======================================================================
 480 //function : SMDS_VolumeTool
 481 //purpose  :
 482 //=======================================================================
 483
 484 SMDS_VolumeTool::~SMDS_VolumeTool()
 485 {
 486   myCurFace.myNodeIndices = NULL;
 487 }
 488
 489 //=======================================================================
 490 //function : SetVolume
 491 //purpose  : Set volume to iterate on
 492 //=======================================================================
 493
 494 bool SMDS_VolumeTool::Set (const SMDS_MeshElement*                  theVolume,
 495                            const bool                               ignoreCentralNodes,
 496                            const std::vector<const SMDS_MeshNode*>* otherNodes)
 497 {
 498   // reset fields
 499   myVolume = 0;
 500   myPolyedre = 0;
 501   myIgnoreCentralNodes = ignoreCentralNodes;
 502
 503   myVolForward = true;
 504   myNbFaces = 0;
 505   myVolumeNodes.clear();
 506   myPolyIndices.clear();
 507   myPolyQuantities.clear();
 508   myPolyFacetOri.clear();
 509   myFwdLinks.clear();
 510
 511   myExternalFaces = false;
 512
 513   myAllFacesNodeIndices_F  = 0;
 514   myAllFacesNodeIndices_RE = 0;
 515   myAllFacesNbNodes        = 0;
 516
 517   myCurFace.myIndex = -1;
 518   myCurFace.myNodeIndices = NULL;
 519   myCurFace.myNodes.clear();
 520
 521   // set volume data
 522   if ( !theVolume || theVolume->GetType() != SMDSAbs_Volume )
 523     return false;
 524
 525   myVolume = theVolume;
 526   myNbFaces = theVolume->NbFaces();
 527   if ( myVolume->IsPoly() )
 528   {
 529     myPolyedre = SMDS_Mesh::DownCast<SMDS_MeshVolume>( myVolume );
 530     myPolyFacetOri.resize( myNbFaces, 0 );
 531   }
 532
 533   // set nodes
 534   myVolumeNodes.resize( myVolume->NbNodes() );
 535   if ( otherNodes )
 536   {
 537     if ( otherNodes->size() != myVolumeNodes.size() )
 538       return ( myVolume = 0 );
 539     for ( size_t i = 0; i < otherNodes->size(); ++i )
 540       if ( ! ( myVolumeNodes[i] = (*otherNodes)[0] ))
 541         return ( myVolume = 0 );
 542   }
 543   else
 544   {
 545     myVolumeNodes.assign( myVolume->begin_nodes(), myVolume->end_nodes() );
 546   }
 547
 548   // check validity
 549   if ( !setFace(0) )
 550     return ( myVolume = 0 );
 551
 552   if ( !myPolyedre )
 553   {
 554     // define volume orientation
 555     XYZ botNormal;
 556     if ( GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z ))
 557     {
 558       const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ];
 559       int topNodeIndex = myVolume->NbCornerNodes() - 1;
 560       while ( !IsLinked( 0, topNodeIndex, /*ignoreMediumNodes=*/true )) --topNodeIndex;
 561       const SMDS_MeshNode* topNode = myVolumeNodes[ topNodeIndex ];
 562       XYZ upDir (topNode->X() - botNode->X(),
 563                  topNode->Y() - botNode->Y(),
 564                  topNode->Z() - botNode->Z() );
 565       myVolForward = ( botNormal.Dot( upDir ) < 0 );
 566     }
 567     if ( !myVolForward )
 568       myCurFace.myIndex = -1; // previous setFace(0) didn't take myVolForward into account
 569   }
 570   return true;
 571 }
 572
 573 //=======================================================================
 574 //function : Inverse
 575 //purpose  : Inverse volume
 576 //=======================================================================
 577
 578 #define SWAP_NODES(nodes,i1,i2)           \
 579 {                                         \
 580   const SMDS_MeshNode* tmp = nodes[ i1 ]; \
 581   nodes[ i1 ] = nodes[ i2 ];              \
 582   nodes[ i2 ] = tmp;                      \
 583 }
 584 void SMDS_VolumeTool::Inverse ()
 585 {
 586   if ( !myVolume ) return;
 587
 588   if (myVolume->IsPoly()) {
 589     MESSAGE("Warning: attempt to inverse polyhedral volume");
 590     return;
 591   }
 592
 593   myVolForward = !myVolForward;
 594   myCurFace.myIndex = -1;
 595
 596   // inverse top and bottom faces
 597   switch ( myVolumeNodes.size() ) {
 598   case 4:
 599     SWAP_NODES( myVolumeNodes, 1, 2 );
 600     break;
 601   case 5:
 602     SWAP_NODES( myVolumeNodes, 1, 3 );
 603     break;
 604   case 6:
 605     SWAP_NODES( myVolumeNodes, 1, 2 );
 606     SWAP_NODES( myVolumeNodes, 4, 5 );
 607     break;
 608   case 8:
 609     SWAP_NODES( myVolumeNodes, 1, 3 );
 610     SWAP_NODES( myVolumeNodes, 5, 7 );
 611     break;
 612   case 12:
 613     SWAP_NODES( myVolumeNodes, 1, 5 );
 614     SWAP_NODES( myVolumeNodes, 2, 4 );
 615     SWAP_NODES( myVolumeNodes, 7, 11 );
 616     SWAP_NODES( myVolumeNodes, 8, 10 );
 617     break;
 618
 619   case 10:
 620     SWAP_NODES( myVolumeNodes, 1, 2 );
 621     SWAP_NODES( myVolumeNodes, 4, 6 );
 622     SWAP_NODES( myVolumeNodes, 8, 9 );
 623     break;
 624   case 13:
 625     SWAP_NODES( myVolumeNodes, 1, 3 );
 626     SWAP_NODES( myVolumeNodes, 5, 8 );
 627     SWAP_NODES( myVolumeNodes, 6, 7 );
 628     SWAP_NODES( myVolumeNodes, 10, 12 );
 629     break;
 630   case 15:
 631     SWAP_NODES( myVolumeNodes, 1, 2 );
 632     SWAP_NODES( myVolumeNodes, 4, 5 );
 633     SWAP_NODES( myVolumeNodes, 6, 8 );
 634     SWAP_NODES( myVolumeNodes, 9, 11 );
 635     SWAP_NODES( myVolumeNodes, 13, 14 );
 636     break;
 637   case 20:
 638     SWAP_NODES( myVolumeNodes, 1, 3 );
 639     SWAP_NODES( myVolumeNodes, 5, 7 );
 640     SWAP_NODES( myVolumeNodes, 8, 11 );
 641     SWAP_NODES( myVolumeNodes, 9, 10 );
 642     SWAP_NODES( myVolumeNodes, 12, 15 );
 643     SWAP_NODES( myVolumeNodes, 13, 14 );
 644     SWAP_NODES( myVolumeNodes, 17, 19 );
 645     break;
 646   case 27:
 647     SWAP_NODES( myVolumeNodes, 1, 3 );
 648     SWAP_NODES( myVolumeNodes, 5, 7 );
 649     SWAP_NODES( myVolumeNodes, 8, 11 );
 650     SWAP_NODES( myVolumeNodes, 9, 10 );
 651     SWAP_NODES( myVolumeNodes, 12, 15 );
 652     SWAP_NODES( myVolumeNodes, 13, 14 );
 653     SWAP_NODES( myVolumeNodes, 17, 19 );
 654     SWAP_NODES( myVolumeNodes, 21, 24 );
 655     SWAP_NODES( myVolumeNodes, 22, 23 );
 656     break;
 657   default:;
 658   }
 659 }
 660
 661 //=======================================================================
 662 //function : GetVolumeType
 663 //purpose  :
 664 //=======================================================================
 665
 666 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetVolumeType() const
 667 {
 668   if ( myPolyedre )
 669     return POLYHEDA;
 670
 671   switch( myVolumeNodes.size() ) {
 672   case 4: return TETRA;
 673   case 5: return PYRAM;
 674   case 6: return PENTA;
 675   case 8: return HEXA;
 676   case 12: return HEX_PRISM;
 677   case 10: return QUAD_TETRA;
 678   case 13: return QUAD_PYRAM;
 679   case 15: return QUAD_PENTA;
 680   case 20: return QUAD_HEXA;
 681   case 27: return QUAD_HEXA;
 682   default: break;
 683   }
 684
 685   return UNKNOWN;
 686 }
 687
 688 //=======================================================================
 689 //function : getTetraVolume
 690 //purpose  :
 691 //=======================================================================
 692
 693 static double getTetraVolume(const SMDS_MeshNode* n1,
 694                              const SMDS_MeshNode* n2,
 695                              const SMDS_MeshNode* n3,
 696                              const SMDS_MeshNode* n4)
 697 {
 698   double p1[3], p2[3], p3[3], p4[3];
 699   n1->GetXYZ( p1 );
 700   n2->GetXYZ( p2 );
 701   n3->GetXYZ( p3 );
 702   n4->GetXYZ( p4 );
 703
 704   double Q1 = -(p1[ 0 ]-p2[ 0 ])*(p3[ 1 ]*p4[ 2 ]-p4[ 1 ]*p3[ 2 ]);
 705   double Q2 =  (p1[ 0 ]-p3[ 0 ])*(p2[ 1 ]*p4[ 2 ]-p4[ 1 ]*p2[ 2 ]);
 706   double R1 = -(p1[ 0 ]-p4[ 0 ])*(p2[ 1 ]*p3[ 2 ]-p3[ 1 ]*p2[ 2 ]);
 707   double R2 = -(p2[ 0 ]-p3[ 0 ])*(p1[ 1 ]*p4[ 2 ]-p4[ 1 ]*p1[ 2 ]);
 708   double S1 =  (p2[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p3[ 2 ]-p3[ 1 ]*p1[ 2 ]);
 709   double S2 = -(p3[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p2[ 2 ]-p2[ 1 ]*p1[ 2 ]);
 710
 711   return (Q1+Q2+R1+R2+S1+S2)/6.0;
 712 }
 713
 714 //=======================================================================
 715 //function : GetSize
 716 //purpose  : Return element volume
 717 //=======================================================================
 718 double SMDS_VolumeTool::GetSize() const
 719 {
 720   double V = 0.;
 721   if ( !myVolume )
 722     return 0.;
 723
 724   if ( myVolume->IsPoly() )
 725   {
 726     if ( !myPolyedre )
 727       return 0.;
 728
 729     SaveFacet savedFacet( myCurFace );
 730
 731     // split a polyhedron into tetrahedrons
 732
 733     bool oriOk = true;
 734     SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
 735     for ( int f = 0; f < NbFaces(); ++f )
 736     {
 737       me->setFace( f );
 738       XYZ area (0,0,0), p1( myCurFace.myNodes[0] );
 739       for ( int n = 0; n < myCurFace.myNbNodes; ++n )
 740       {
 741         XYZ p2( myCurFace.myNodes[ n+1 ]);
 742         area = area + p1.Crossed( p2 );
 743         p1 = p2;
 744       }
 745       V += p1.Dot( area );
 746       oriOk = oriOk && IsFaceExternal( f );
 747     }
 748     V /= 6;
 749     if ( !oriOk && V > 0 )
 750       V *= -1;
 751   }
 752   else
 753   {
 754     const static int ind[] = {
 755       0, 1, 3, 6, 11, 23, 31, 44, 58, 78 };
 756     const static int vtab[][4] = { // decomposition into tetra in the order of enum VolumeType
 757       // tetrahedron
 758       { 0, 1, 2, 3 },
 759       // pyramid
 760       { 0, 1, 3, 4 },
 761       { 1, 2, 3, 4 },
 762       // pentahedron
 763       { 0, 1, 2, 3 },
 764       { 1, 5, 3, 4 },
 765       { 1, 5, 2, 3 },
 766       // hexahedron
 767       { 1, 4, 3, 0 },
 768       { 4, 1, 6, 5 },
 769       { 1, 3, 6, 2 },
 770       { 4, 6, 3, 7 },
 771       { 1, 4, 6, 3 },
 772       // hexagonal prism
 773       { 0, 1, 2, 7 },
 774       { 0, 7, 8, 6 },
 775       { 2, 7, 8, 0 },
 776
 777       { 0, 3, 4, 9 },
 778       { 0, 9, 10, 6 },
 779       { 4, 9, 10, 0 },
 780
 781       { 0, 3, 4, 9 },
 782       { 0, 9, 10, 6 },
 783       { 4, 9, 10, 0 },
 784
 785       { 0, 4, 5, 10 },
 786       { 0, 10, 11, 6 },
 787       { 5, 10, 11, 0 },
 788
 789       // quadratic tetrahedron
 790       { 0, 4, 6, 7 },
 791       { 1, 5, 4, 8 },
 792       { 2, 6, 5, 9 },
 793       { 7, 8, 9, 3 },
 794       { 4, 6, 7, 9 },
 795       { 4, 5, 6, 9 },
 796       { 4, 7, 8, 9 },
 797       { 4, 5, 9, 8 },
 798
 799       // quadratic pyramid
 800       { 0, 5, 8, 9 },
 801       { 1, 5,10, 6 },
 802       { 2, 6,11, 7 },
 803       { 3, 7,12, 8 },
 804       { 4, 9,11,10 },
 805       { 4, 9,12,11 },
 806       { 10, 5, 9, 8 },
 807       { 10, 8, 9,12 },
 808       { 10, 8,12, 7 },
 809       { 10, 7,12,11 },
 810       { 10, 7,11, 6 },
 811       { 10, 5, 8, 6 },
 812       { 10, 6, 8, 7 },
 813
 814       // quadratic pentahedron
 815       { 12, 0, 8, 6 },
 816       { 12, 8, 7, 6 },
 817       { 12, 8, 2, 7 },
 818       { 12, 6, 7, 1 },
 819       { 12, 1, 7,13 },
 820       { 12, 7, 2,13 },
 821       { 12, 2,14,13 },
 822
 823       { 12, 3, 9,11 },
 824       { 12,11, 9,10 },
 825       { 12,11,10, 5 },
 826       { 12, 9, 4,10 },
 827       { 12,14, 5,10 },
 828       { 12,14,10, 4 },
 829       { 12,14, 4,13 },
 830
 831       // quadratic hexahedron
 832       { 16, 0,11, 8 },
 833       { 16,11, 9, 8 },
 834       { 16, 8, 9, 1 },
 835       { 16,11, 3,10 },
 836       { 16,11,10, 9 },
 837       { 16,10, 2, 9 },
 838       { 16, 3,19, 2 },
 839       { 16, 2,19,18 },
 840       { 16, 2,18,17 },
 841       { 16, 2,17, 1 },
 842
 843       { 16, 4,12,15 },
 844       { 16,12, 5,13 },
 845       { 16,12,13,15 },
 846       { 16,13, 6,14 },
 847       { 16,13,14,15 },
 848       { 16,14, 7,15 },
 849       { 16, 6, 5,17 },
 850       { 16,18, 6,17 },
 851       { 16,18, 7, 6 },
 852       { 16,18,19, 7 },
 853
 854     };
 855
 856     int type = GetVolumeType();
 857     int n1 = ind[type];
 858     int n2 = ind[type+1];
 859
 860     for (int i = n1; i <  n2; i++) {
 861       V -= getTetraVolume( myVolumeNodes[ vtab[i][0] ],
 862                            myVolumeNodes[ vtab[i][1] ],
 863                            myVolumeNodes[ vtab[i][2] ],
 864                            myVolumeNodes[ vtab[i][3] ]);
 865     }
 866   }
 867   return V;
 868 }
 869
 870
 871 //=======================================================================
 872 //function : getTetraScaledJacobian
 873 //purpose  : Given the smesh nodes in the canonical order of the tetrahedron, return the scaled jacobian
 874 //=======================================================================
 875 static double getTetraScaledJacobian(const SMDS_MeshNode* n0,
 876                                      const SMDS_MeshNode* n1,
 877                                      const SMDS_MeshNode* n2,
 878                                      const SMDS_MeshNode* n3)
 879 {
 880   const double sqrt = std::sqrt(2.0);
 881   // Get the coordinates
 882   XYZ p0( n0 );
 883   XYZ p1( n1 );
 884   XYZ p2( n2 );
 885   XYZ p3( n3 );
 886   // Define the edges connecting the nodes
 887   XYZ L0 = p1-p0;
 888   XYZ L1 = p2-p1;
 889   XYZ L2 = p2-p0; // invert the definition of doc to get the proper orientation of the crossed product
 890   XYZ L3 = p3-p0;
 891   XYZ L4 = p3-p1;
 892   XYZ L5 = p3-p2;
 893   double Jacobian = L2.Crossed( L0 ).Dot( L3 );
 894   double norm0 = L0.Magnitude();
 895   double norm1 = L1.Magnitude();
 896   double norm2 = L2.Magnitude();
 897   double norm3 = L3.Magnitude();
 898   double norm4 = L4.Magnitude();
 899   double norm5 = L5.Magnitude();
 900
 901   std::array<double, 5> norms{};
 902   norms[0] = Jacobian;
 903   norms[1] = norm3*norm4*norm5;
 904   norms[2] = norm1*norm2*norm5;
 905   norms[3] = norm0*norm1*norm4;
 906   norms[4] = norm0*norm2*norm3;
 907
 908   auto findMaxNorm = std::max_element(norms.begin(), norms.end());
 909   double maxNorm = *findMaxNorm;
 910
 911   if ( std::fabs( maxNorm ) < std::numeric_limits<double>::min() )
 912     maxNorm = std::numeric_limits<double>::max();
 913
 914   return Jacobian * sqrt / maxNorm;
 915 }
 916
 917 //=======================================================================
 918 //function : getPyramidScaledJacobian
 919 //purpose  : Given the pyramid, compute the scaled jacobian of the four tetrahedrons and return the minimun value.
 920 //=======================================================================
 921 static double getPyramidScaledJacobian(const SMDS_MeshNode* n0,
 922                                         const SMDS_MeshNode* n1,
 923                                         const SMDS_MeshNode* n2,
 924                                         const SMDS_MeshNode* n3,
 925                                         const SMDS_MeshNode* n4)
 926 {
 927   const double sqrt = std::sqrt(2.0);
 928   std::array<double, 4> tetScaledJacobian{};
 929   tetScaledJacobian[0] = getTetraScaledJacobian(n0, n1, n3, n4);
 930   tetScaledJacobian[1] = getTetraScaledJacobian(n1, n2, n0, n4);
 931   tetScaledJacobian[2] = getTetraScaledJacobian(n2, n3, n1, n4);
 932   tetScaledJacobian[3] = getTetraScaledJacobian(n3, n0, n2, n4);
 933
 934   auto minEntry = std::min_element(tetScaledJacobian.begin(), tetScaledJacobian.end());
 935
 936   double scaledJacobian = (*minEntry) * 2.0/sqrt;
 937   return scaledJacobian < 1.0 ? scaledJacobian : 1.0 - (scaledJacobian - 1.0);
 938 }
 939
 940
 941
 942 //=======================================================================
 943 //function : getHexaScaledJacobian
 944 //purpose  : Evaluate the scaled jacobian on the eight vertices of the hexahedron and return the minimal registered value
 945 //remark   : Follow the reference numeration described at the top of the class.
 946 //=======================================================================
 947 static double getHexaScaledJacobian(const SMDS_MeshNode* n0,
 948                                     const SMDS_MeshNode* n1,
 949                                     const SMDS_MeshNode* n2,
 950                                     const SMDS_MeshNode* n3,
 951                                     const SMDS_MeshNode* n4,
 952                                     const SMDS_MeshNode* n5,
 953                                     const SMDS_MeshNode* n6,
 954                                     const SMDS_MeshNode* n7)
 955 {
 956   // Scaled jacobian is an scalar quantity measuring the deviation of the geometry from the perfect geometry
 957   // Get the coordinates
 958   XYZ p0( n0 );
 959   XYZ p1( n1 );
 960   XYZ p2( n2 );
 961   XYZ p3( n3 );
 962   XYZ p4( n4 );
 963   XYZ p5( n5 );
 964   XYZ p6( n6 );
 965   XYZ p7( n7 );
 966
 967   // Define the edges connecting the nodes
 968   XYZ L0  = (p1-p0).Normalize();
 969   XYZ L1  = (p2-p1).Normalize();
 970   XYZ L2  = (p3-p2).Normalize();
 971   XYZ L3  = (p3-p0).Normalize();
 972   XYZ L4  = (p4-p0).Normalize();
 973   XYZ L5  = (p5-p1).Normalize();
 974   XYZ L6  = (p6-p2).Normalize();
 975   XYZ L7  = (p7-p3).Normalize();
 976   XYZ L8  = (p5-p4).Normalize();
 977   XYZ L9  = (p6-p5).Normalize();
 978   XYZ L10 = (p7-p6).Normalize();
 979   XYZ L11 = (p7-p4).Normalize();
 980   XYZ X0  = (p1-p0+p2-p3+p6-p7+p5-p4).Normalize();
 981   XYZ X1  = (p3-p0+p2-p1+p7-p4+p6-p5).Normalize();
 982   XYZ X2  = (p4-p0+p7-p3+p5-p1+p6-p2).Normalize();
 983
 984   std::array<double, 9> scaledJacobian{};
 985   //Scaled jacobian of nodes following their numeration
 986   scaledJacobian[0] =  L4.Crossed( L3).Dot( L0 );   // For L0
 987   scaledJacobian[1] =  L5.Crossed(-L0).Dot( L1 );   // For L1
 988   scaledJacobian[2] =  L6.Crossed(-L1).Dot( L2 );   // For L2
 989   scaledJacobian[3] =  L7.Crossed(-L2).Dot(-L3 );   // For L3
 990   scaledJacobian[4] = -L4.Crossed( L8).Dot( L11 );  // For L11
 991   scaledJacobian[5] = -L5.Crossed( L9).Dot(-L8 );   // For L8
 992   scaledJacobian[6] = -L6.Crossed(L10).Dot(-L9 );   // For L9
 993   scaledJacobian[7] = -L7.Crossed(-L11).Dot(-L10 ); // For L10
 994   scaledJacobian[8] =  X2.Crossed( X1).Dot( X0 );   // For principal axes
 995
 996   auto minScaledJacobian = std::min_element(scaledJacobian.begin(), scaledJacobian.end());
 997   return *minScaledJacobian;
 998 }
 999
1000
1001 //=======================================================================
1002 //function : getTetraNormalizedJacobian
1003 //purpose  : Return the jacobian of the tetrahedron based on normalized vectors
1004 //=======================================================================
1005 static double getTetraNormalizedJacobian(const SMDS_MeshNode* n0,
1006                                           const SMDS_MeshNode* n1,
1007                                           const SMDS_MeshNode* n2,
1008                                           const SMDS_MeshNode* n3)
1009 {
1010   const double sqrt = std::sqrt(2.0);
1011   // Get the coordinates
1012   XYZ p0( n0 );
1013   XYZ p1( n1 );
1014   XYZ p2( n2 );
1015   XYZ p3( n3 );
1016   // Define the normalized edges connecting the nodes
1017   XYZ L0 = (p1-p0).Normalize();
1018   XYZ L2 = (p2-p0).Normalize(); // invert the definition of doc to get the proper orientation of the crossed product
1019   XYZ L3 = (p3-p0).Normalize();
1020   return L2.Crossed( L0 ).Dot( L3 );
1021 }
1022
1023 //=======================================================================
1024 //function : getPentaScaledJacobian
1025 //purpose  : Evaluate the scaled jacobian on the pentahedron based on decomposed tetrahedrons
1026 //=======================================================================
1027 /*
1028 //            + N1
1029 //           /|\
1030 //          / | \
1031 //         /  |  \
1032 //        /   |   \
1033 //    N0 +---------+ N2
1034 //       |    |    |               NUMERATION RERENCE FOLLOWING POSSITIVE RIGHT HAND RULE
1035 //       |    + N4 |
1036 //       |   / \   |               PENTAHEDRON
1037 //       |  /   \  |
1038 //       | /     \ |
1039 //       |/       \|
1040 //    N3 +---------+ N5
1041 //
1042 //          N1
1043 //          +
1044 //          |\
1045 //         /| \
1046 //        / |  \
1047 //    N0 +--|---+ N2               TETRAHEDRON ASSOCIATED TO N0
1048 //       \  |   /                  Numeration passed to getTetraScaledJacobian
1049 //        \ |  /                   N0=N0; N1=N2; N2=N3; N3=N1
1050 //         \| /
1051 //          |/
1052 //          +
1053 //          N3
1054 //
1055 //           N1
1056 //           +
1057 //          /|\
1058 //         / | \
1059 //        /  |  \
1060 //    N2 +---|---+ N5             TETRAHEDRON ASSOCIATED TO N2
1061 //       \   |  /                 Numeration passed to getTetraScaledJacobian
1062 //        \  | /                  N0=N2; N1=N5; N2=N0; N3=N1
1063 //         \ |/
1064 //          \|
1065 //           +
1066 //           N0
1067 //
1068 //           N4
1069 //           +
1070 //          /|\
1071 //         / | \
1072 //        /  |  \
1073 //    N3 +---|---+ N0             TETRAHEDRON ASSOCIATED TO N3
1074 //       \   |   /                Numeration passed to getTetraScaledJacobian
1075 //        \  |  /                 N0=N3; N1=N0; N2=N5; N3=N4
1076 //         \ | /
1077 //          \|/
1078 //           +
1079 //           N5
1080 //
1081 //           N3
1082 //           +
1083 //          /|\
1084 //         / | \
1085 //        /  |  \
1086 //    N1 +---|---+ N2             TETRAHEDRON ASSOCIATED TO N1
1087 //       \   |   /                Numeration passed to getTetraScaledJacobian
1088 //        \  |  /                 N0=N1; N1=N2; N2=N0; N3=N3
1089 //         \ | /
1090 //          \|/
1091 //           +
1092 //           N0
1093 //
1094 //          ...
1095 */
1096 static double getPentaScaledJacobian(const SMDS_MeshNode* n0,
1097                                      const SMDS_MeshNode* n1,
1098                                      const SMDS_MeshNode* n2,
1099                                      const SMDS_MeshNode* n3,
1100                                      const SMDS_MeshNode* n4,
1101                                      const SMDS_MeshNode* n5)
1102 {
1103   std::array<double, 6> scaledJacobianOfReferenceTetra{};
1104   scaledJacobianOfReferenceTetra[0] = getTetraNormalizedJacobian(n0, n2, n3, n1);  // For n0
1105   scaledJacobianOfReferenceTetra[1] = getTetraNormalizedJacobian(n2, n5, n0, n1);  // For n2
1106   scaledJacobianOfReferenceTetra[2] = getTetraNormalizedJacobian(n3, n0, n5, n4);  // For n3
1107   scaledJacobianOfReferenceTetra[3] = getTetraNormalizedJacobian(n5, n3, n2, n4);  // For n5
1108   scaledJacobianOfReferenceTetra[4] = getTetraNormalizedJacobian(n1, n2, n0, n3);  // For n1
1109   scaledJacobianOfReferenceTetra[5] = getTetraNormalizedJacobian(n4, n3, n5, n2);  // For n4
1110
1111   auto minScaledJacobian = std::min_element(scaledJacobianOfReferenceTetra.begin(), scaledJacobianOfReferenceTetra.end());
1112   double minScalJac = (*minScaledJacobian)* 2.0 / std::sqrt(3.0);
1113
1114   if (minScalJac > 0)
1115     return std::min(minScalJac, std::numeric_limits<double>::max());
1116
1117   return std::max(minScalJac, -std::numeric_limits<double>::max());
1118 }
1119
1120 //=======================================================================
1121 //function : getHexaPrismScaledJacobian
1122 //purpose  : Evaluate the scaled jacobian on the hexaprism by decomposing the goemetry into three 1hexa + 2 pentahedrons
1123 //=======================================================================
1124 static double getHexaPrismScaledJacobian(const SMDS_MeshNode* n0,
1125                                           const SMDS_MeshNode* n1,
1126                                           const SMDS_MeshNode* n2,
1127                                           const SMDS_MeshNode* n3,
1128                                           const SMDS_MeshNode* n4,
1129                                           const SMDS_MeshNode* n5,
1130                                           const SMDS_MeshNode* n6,
1131                                           const SMDS_MeshNode* n7,
1132                                           const SMDS_MeshNode* n8,
1133                                           const SMDS_MeshNode* n9,
1134                                           const SMDS_MeshNode* n10,
1135                                           const SMDS_MeshNode* n11)
1136 {
1137   // The Pentahedron from the left
1138   // n0=n0; n1=n1; n2=n2; n3=n6; n4=n7, n5=n8;
1139   double scaledJacobianPentleft = getPentaScaledJacobian( n0, n1, n2, n6, n7, n8 );
1140   // The core Hexahedron
1141   // n0=n0; n1=n2, n2=n3; n3=n5; n4=n6; n5=n8; n6=n9; n7=n11
1142   double scaledJacobianHexa     = getHexaScaledJacobian( n0, n2, n3, n5, n6, n8, n9, n11 );
1143   // The Pentahedron from the right
1144   // n0=n5; n1=n4; n2=n3; n3=n11; n4=n10; n5=n9
1145   double scaledJacobianPentright = getPentaScaledJacobian( n5, n4, n3, n11, n10, n9 );
1146
1147   return std::min( scaledJacobianHexa, std::min( scaledJacobianPentleft, scaledJacobianPentright ) );
1148
1149 }
1150
1151 //=======================================================================
1152 //function : GetScaledJacobian
1153 //purpose  : Return element Scaled Jacobian using the generic definition given
1154 //            in https://gitlab.kitware.com/third-party/verdict/-/blob/master/SAND2007-2853p.pdf
1155 //=======================================================================
1156
1157 double SMDS_VolumeTool::GetScaledJacobian() const
1158 {
1159
1160   // For Tetra, call directly the getTetraScaledJacobian
1161   double scaledJacobian = 0.;
1162
1163   VolumeType type = GetVolumeType();
1164   switch (type)
1165   {
1166   case TETRA:
1167   case QUAD_TETRA:
1168     scaledJacobian = getTetraScaledJacobian( myVolumeNodes[0], myVolumeNodes[1], myVolumeNodes[2], myVolumeNodes[3] );
1169     break;
1170   case HEXA:
1171   case QUAD_HEXA:
1172     scaledJacobian = getHexaScaledJacobian( myVolumeNodes[0], myVolumeNodes[1], myVolumeNodes[2], myVolumeNodes[3],
1173                                             myVolumeNodes[4], myVolumeNodes[5], myVolumeNodes[6], myVolumeNodes[7] );
1174     break;
1175   case PYRAM:
1176   case QUAD_PYRAM:
1177     scaledJacobian = getPyramidScaledJacobian( myVolumeNodes[0], myVolumeNodes[1], myVolumeNodes[2], myVolumeNodes[3], myVolumeNodes[4] );
1178     break;
1179   case PENTA:
1180   case QUAD_PENTA:
1181     scaledJacobian = getPentaScaledJacobian( myVolumeNodes[0], myVolumeNodes[1],
1182                                              myVolumeNodes[2], myVolumeNodes[3],
1183                                              myVolumeNodes[4], myVolumeNodes[5] );
1184     break;
1185   case HEX_PRISM:
1186     scaledJacobian = getHexaPrismScaledJacobian( myVolumeNodes[0], myVolumeNodes[1], myVolumeNodes[2], myVolumeNodes[3],
1187                                                  myVolumeNodes[4], myVolumeNodes[5], myVolumeNodes[6], myVolumeNodes[7],
1188                                                  myVolumeNodes[8], myVolumeNodes[9], myVolumeNodes[10], myVolumeNodes[11]);
1189     break;
1190   default:
1191     break;
1192   }
1193
1194   return scaledJacobian;
1195 }
1196
1197
1198 //=======================================================================
1199 //function : GetBaryCenter
1200 //purpose  :
1201 //=======================================================================
1202
1203 bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
1204 {
1205   X = Y = Z = 0.;
1206   if ( !myVolume )
1207     return false;
1208
1209   for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
1210     X += myVolumeNodes[ i ]->X();
1211     Y += myVolumeNodes[ i ]->Y();
1212     Z += myVolumeNodes[ i ]->Z();
1213   }
1214   X /= myVolumeNodes.size();
1215   Y /= myVolumeNodes.size();
1216   Z /= myVolumeNodes.size();
1217
1218   return true;
1219 }
1220
1221 //================================================================================
1222 /*!
1223  * \brief Classify a point
1224  *  \param tol - thickness of faces
1225  */
1226 //================================================================================
1227
1228 bool SMDS_VolumeTool::IsOut(double X, double Y, double Z, double tol) const
1229 {
1230   // LIMITATION: for convex volumes only
1231   XYZ p( X,Y,Z );
1232   for ( int iF = 0; iF < myNbFaces; ++iF )
1233   {
1234     XYZ faceNormal;
1235     if ( !GetFaceNormal( iF, faceNormal.x, faceNormal.y, faceNormal.z ))
1236       continue;
1237     if ( !IsFaceExternal( iF ))
1238       faceNormal = XYZ() - faceNormal; // reverse
1239
1240     XYZ face2p( p - XYZ( myCurFace.myNodes[0] ));
1241     if ( face2p.Dot( faceNormal ) > tol )
1242       return true;
1243   }
1244   return false;
1245 }
1246
1247 //=======================================================================
1248 //function : SetExternalNormal
1249 //purpose  : Node order will be so that faces normals are external
1250 //=======================================================================
1251
1252 void SMDS_VolumeTool::SetExternalNormal ()
1253 {
1254   myExternalFaces = true;
1255   myCurFace.myIndex = -1;
1256 }
1257
1258 //=======================================================================
1259 //function : NbFaceNodes
1260 //purpose  : Return number of nodes in the array of face nodes
1261 //=======================================================================
1262
1263 int SMDS_VolumeTool::NbFaceNodes( int faceIndex ) const
1264 {
1265   if ( !setFace( faceIndex ))
1266     return 0;
1267   return myCurFace.myNbNodes;
1268 }
1269
1270 //=======================================================================
1271 //function : GetFaceNodes
1272 //purpose  : Return pointer to the array of face nodes.
1273 //           To comfort link iteration, the array
1274 //           length == NbFaceNodes( faceIndex ) + 1 and
1275 //           the last node == the first one.
1276 //=======================================================================
1277
1278 const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex ) const
1279 {
1280   if ( !setFace( faceIndex ))
1281     return 0;
1282   return &myCurFace.myNodes[0];
1283 }
1284
1285 //=======================================================================
1286 //function : GetFaceNodesIndices
1287 //purpose  : Return pointer to the array of face nodes indices
1288 //           To comfort link iteration, the array
1289 //           length == NbFaceNodes( faceIndex ) + 1 and
1290 //           the last node index == the first one.
1291 //=======================================================================
1292
1293 const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex ) const
1294 {
1295   if ( !setFace( faceIndex ))
1296     return 0;
1297
1298   return myCurFace.myNodeIndices;
1299 }
1300
1301 //=======================================================================
1302 //function : GetFaceNodes
1303 //purpose  : Return a set of face nodes.
1304 //=======================================================================
1305
1306 bool SMDS_VolumeTool::GetFaceNodes (int                             faceIndex,
1307                                     std::set<const SMDS_MeshNode*>& theFaceNodes ) const
1308 {
1309   if ( !setFace( faceIndex ))
1310     return false;
1311
1312   theFaceNodes.clear();
1313   theFaceNodes.insert( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
1314
1315   return true;
1316 }
1317
1318 namespace
1319 {
1320   struct NLink : public std::pair<smIdType,smIdType>
1321   {
1322     int myOri;
1323     NLink(const SMDS_MeshNode* n1=0, const SMDS_MeshNode* n2=0, int ori=1 )
1324     {
1325       if ( n1 )
1326       {
1327         if (( myOri = ( n1->GetID() < n2->GetID() )))
1328         {
1329           first  = n1->GetID();
1330           second = n2->GetID();
1331         }
1332         else
1333         {
1334           myOri  = -1;
1335           first  = n2->GetID();
1336           second = n1->GetID();
1337         }
1338         myOri *= ori;
1339       }
1340       else
1341       {
1342         myOri = first = second = 0;
1343       }
1344     }
1345     //int Node1() const { return myOri == -1 ? second : first; }
1346
1347     //bool IsSameOri( const std::pair<int,int>& link ) const { return link.first == Node1(); }
1348   };
1349 }
1350
1351 //=======================================================================
1352 //function : IsFaceExternal
1353 //purpose  : Check normal orientation of a given face
1354 //=======================================================================
1355
1356 bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) const
1357 {
1358   if ( myExternalFaces || !myVolume )
1359     return true;
1360
1361   if ( !myPolyedre ) // all classical volumes have external facet normals
1362     return true;
1363
1364   SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
1365
1366   if ( myPolyFacetOri[ faceIndex ])
1367     return myPolyFacetOri[ faceIndex ] > 0;
1368
1369   int ori = 0; // -1-in, +1-out, 0-undef
1370   double minProj, maxProj;
1371   if ( projectNodesToNormal( faceIndex, minProj, maxProj ))
1372   {
1373     // all nodes are on the same side of the facet
1374     ori = ( minProj < 0 ? +1 : -1 );
1375     me->myPolyFacetOri[ faceIndex ] = ori;
1376
1377     if ( !myFwdLinks.empty() ) // concave polyhedron; collect oriented links
1378       for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1379       {
1380         NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
1381         me->myFwdLinks.insert( make_pair( link, link.myOri ));
1382       }
1383     return ori > 0;
1384   }
1385
1386   SaveFacet savedFacet( myCurFace );
1387
1388   // concave polyhedron
1389
1390   if ( myFwdLinks.empty() ) // get links of the least ambiguously oriented facet
1391   {
1392     for ( size_t i = 0; i < myPolyFacetOri.size() && !ori; ++i )
1393       ori = myPolyFacetOri[ i ];
1394
1395     if ( !ori ) // none facet is oriented yet
1396     {
1397       // find the least ambiguously oriented facet
1398       int faceMostConvex = -1;
1399       std::map< double, int > convexity2face;
1400       for ( size_t iF = 0; iF < myPolyFacetOri.size() && faceMostConvex < 0; ++iF )
1401       {
1402         if ( projectNodesToNormal( iF, minProj, maxProj ))
1403         {
1404           // all nodes are on the same side of the facet
1405           me->myPolyFacetOri[ iF ] = ( minProj < 0 ? +1 : -1 );
1406           faceMostConvex = iF;
1407         }
1408         else
1409         {
1410           ori = ( -minProj < maxProj ? -1 : +1 );
1411           double convexity = std::min( -minProj, maxProj ) / std::max( -minProj, maxProj );
1412           convexity2face.insert( std::make_pair( convexity, iF * ori ));
1413         }
1414       }
1415       if ( faceMostConvex < 0 ) // none facet has nodes on the same side
1416       {
1417         // use the least ambiguous facet
1418         faceMostConvex = convexity2face.begin()->second;
1419         ori = ( faceMostConvex < 0 ? -1 : +1 );
1420         faceMostConvex = std::abs( faceMostConvex );
1421         me->myPolyFacetOri[ faceMostConvex ] = ori;
1422       }
1423     }
1424     // collect links of the oriented facets in myFwdLinks
1425     for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
1426     {
1427       ori = myPolyFacetOri[ iF ];
1428       if ( !ori ) continue;
1429       setFace( iF );
1430       for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1431       {
1432         NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
1433         me->myFwdLinks.insert( make_pair( link, link.myOri ));
1434       }
1435     }
1436   }
1437
1438   // compare orientation of links of the facet with myFwdLinks
1439   ori = 0;
1440   setFace( faceIndex );
1441   std::vector< NLink > links( myCurFace.myNbNodes ), links2;
1442   for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
1443   {
1444     NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
1445     std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
1446     if ( l2o != myFwdLinks.end() )
1447       ori = link.myOri * l2o->second * -1;
1448     links[ i ] = link;
1449   }
1450   while ( !ori ) // the facet has no common links with already oriented facets
1451   {
1452     // orient and collect links of other non-oriented facets
1453     for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
1454     {
1455       if ( myPolyFacetOri[ iF ] ) continue; // already oriented
1456       setFace( iF );
1457       links2.clear();
1458       ori = 0;
1459       for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
1460       {
1461         NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
1462         std::map<Link, int>::const_iterator l2o = myFwdLinks.find( link );
1463         if ( l2o != myFwdLinks.end() )
1464           ori = link.myOri * l2o->second * -1;
1465         links2.push_back( link );
1466       }
1467       if ( ori ) // one more facet oriented
1468       {
1469         me->myPolyFacetOri[ iF ] = ori;
1470         for ( size_t i = 0; i < links2.size(); ++i )
1471           me->myFwdLinks.insert( make_pair( links2[i], links2[i].myOri * ori ));
1472         break;
1473       }
1474     }
1475     if ( !ori )
1476       return false; // error in algorithm: infinite loop
1477
1478     // try to orient the facet again
1479     ori = 0;
1480     for ( size_t i = 0; i < links.size() && !ori; ++i )
1481     {
1482       std::map<Link, int>::const_iterator l2o = myFwdLinks.find( links[i] );
1483       if ( l2o != myFwdLinks.end() )
1484         ori = links[i].myOri * l2o->second * -1;
1485     }
1486     me->myPolyFacetOri[ faceIndex ] = ori;
1487   }
1488
1489   return ori > 0;
1490 }
1491
1492 //=======================================================================
1493 //function : projectNodesToNormal
1494 //purpose  : compute min and max projections of all nodes to normal of a facet.
1495 //=======================================================================
1496
1497 bool SMDS_VolumeTool::projectNodesToNormal( int     faceIndex,
1498                                             double& minProj,
1499                                             double& maxProj,
1500                                             double* normalXYZ ) const
1501 {
1502   minProj = std::numeric_limits<double>::max();
1503   maxProj = std::numeric_limits<double>::min();
1504
1505   XYZ normal;
1506   if ( !GetFaceNormal( faceIndex, normal.x, normal.y, normal.z ))
1507     return false;
1508   if ( normalXYZ )
1509     memcpy( normalXYZ, normal.data(), 3*sizeof(double));
1510
1511   XYZ p0 ( myCurFace.myNodes[0] );
1512   for ( size_t i = 0; i < myVolumeNodes.size(); ++i )
1513   {
1514     if ( std::find( myCurFace.myNodes.begin() + 1,
1515                     myCurFace.myNodes.end(),
1516                     myVolumeNodes[ i ] ) != myCurFace.myNodes.end() )
1517       continue; // node of the faceIndex-th facet
1518
1519     double proj = normal.Dot( XYZ( myVolumeNodes[ i ]) - p0 );
1520     if ( proj < minProj ) minProj = proj;
1521     if ( proj > maxProj ) maxProj = proj;
1522   }
1523   const double tol = 1e-7;
1524   minProj += tol;
1525   maxProj -= tol;
1526   bool diffSize = ( minProj * maxProj < 0 );
1527   // if ( diffSize )
1528   // {
1529   //   minProj = -minProj;
1530   // }
1531   // else if ( minProj < 0 )
1532   // {
1533   //   minProj = -minProj;
1534   //   maxProj = -maxProj;
1535   // }
1536
1537   return !diffSize; // ? 0 : (minProj >= 0);
1538 }
1539
1540 //=======================================================================
1541 //function : GetFaceNormal
1542 //purpose  : Return a normal to a face
1543 //=======================================================================
1544
1545 bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z) const
1546 {
1547   if ( !setFace( faceIndex ))
1548     return false;
1549
1550   const int iQuad = ( !myPolyedre && myCurFace.myNbNodes > 6 ) ? 2 : 1;
1551   XYZ p1 ( myCurFace.myNodes[0*iQuad] );
1552   XYZ p2 ( myCurFace.myNodes[1*iQuad] );
1553   XYZ p3 ( myCurFace.myNodes[2*iQuad] );
1554   XYZ aVec12( p2 - p1 );
1555   XYZ aVec13( p3 - p1 );
1556   XYZ cross = aVec12.Crossed( aVec13 );
1557
1558   for ( int i = 3*iQuad; i < myCurFace.myNbNodes; i += iQuad )
1559   {
1560     XYZ p4 ( myCurFace.myNodes[i] );
1561     XYZ aVec14( p4 - p1 );
1562     XYZ cross2 = aVec13.Crossed( aVec14 );
1563     cross = cross + cross2;
1564     aVec13 = aVec14;
1565   }
1566
1567   double size = cross.Magnitude();
1568   if ( size <= std::numeric_limits<double>::min() )
1569     return false;
1570
1571   X = cross.x / size;
1572   Y = cross.y / size;
1573   Z = cross.z / size;
1574
1575   return true;
1576 }
1577
1578 //================================================================================
1579 /*!
1580  * \brief Return barycenter of a face
1581  */
1582 //================================================================================
1583
1584 bool SMDS_VolumeTool::GetFaceBaryCenter (int faceIndex, double & X, double & Y, double & Z) const
1585 {
1586   if ( !setFace( faceIndex ))
1587     return false;
1588
1589   X = Y = Z = 0.0;
1590   for ( int i = 0; i < myCurFace.myNbNodes; ++i )
1591   {
1592     X += myCurFace.myNodes[i]->X() / myCurFace.myNbNodes;
1593     Y += myCurFace.myNodes[i]->Y() / myCurFace.myNbNodes;
1594     Z += myCurFace.myNodes[i]->Z() / myCurFace.myNbNodes;
1595   }
1596   return true;
1597 }
1598
1599 //=======================================================================
1600 //function : GetFaceArea
1601 //purpose  : Return face area
1602 //=======================================================================
1603
1604 double SMDS_VolumeTool::GetFaceArea( int faceIndex ) const
1605 {
1606   double area = 0;
1607   if ( !setFace( faceIndex ))
1608     return area;
1609
1610   XYZ p1 ( myCurFace.myNodes[0] );
1611   XYZ p2 ( myCurFace.myNodes[1] );
1612   XYZ p3 ( myCurFace.myNodes[2] );
1613   XYZ aVec12( p2 - p1 );
1614   XYZ aVec13( p3 - p1 );
1615   area += aVec12.Crossed( aVec13 ).Magnitude();
1616
1617   if (myVolume->IsPoly())
1618   {
1619     for ( int i = 3; i < myCurFace.myNbNodes; ++i )
1620     {
1621       XYZ pI ( myCurFace.myNodes[i] );
1622       XYZ aVecI( pI - p1 );
1623       area += aVec13.Crossed( aVecI ).Magnitude();
1624       aVec13 = aVecI;
1625     }
1626   }
1627   else
1628   {
1629     if ( myCurFace.myNbNodes == 4 ) {
1630       XYZ p4 ( myCurFace.myNodes[3] );
1631       XYZ aVec14( p4 - p1 );
1632       area += aVec14.Crossed( aVec13 ).Magnitude();
1633     }
1634   }
1635   return area / 2;
1636 }
1637
1638 //================================================================================
1639 /*!
1640  * \brief Return index of the node located at face center of a quadratic element like HEX27
1641  */
1642 //================================================================================
1643
1644 int SMDS_VolumeTool::GetCenterNodeIndex( int faceIndex ) const
1645 {
1646   if ( myAllFacesNbNodes && myVolumeNodes.size() == 27 ) // classic element with 27 nodes
1647   {
1648     switch ( faceIndex ) {
1649     case 0: return 20;
1650     case 1: return 25;
1651     default:
1652       return faceIndex + 19;
1653     }
1654   }
1655   return -1;
1656 }
1657
1658 //=======================================================================
1659 //function : GetOppFaceIndex
1660 //purpose  : Return index of the opposite face if it exists, else -1.
1661 //=======================================================================
1662
1663 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
1664 {
1665   int ind = -1;
1666   if (myPolyedre) {
1667     MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume");
1668     return ind;
1669   }
1670
1671   const int nbHoriFaces = 2;
1672
1673   if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
1674     switch ( myVolumeNodes.size() ) {
1675     case 6:
1676     case 15:
1677       if ( faceIndex == 0 || faceIndex == 1 )
1678         ind = 1 - faceIndex;
1679       break;
1680     case 8:
1681     case 12:
1682       if ( faceIndex <= 1 ) // top or bottom
1683         ind = 1 - faceIndex;
1684       else {
1685         const int nbSideFaces = myAllFacesNbNodes[0];
1686         ind = ( faceIndex - nbHoriFaces + nbSideFaces/2 ) % nbSideFaces + nbHoriFaces;
1687       }
1688       break;
1689     case 20:
1690     case 27:
1691       ind = GetOppFaceIndexOfHex( faceIndex );
1692       break;
1693     default:;
1694     }
1695   }
1696   return ind;
1697 }
1698
1699 //=======================================================================
1700 //function : GetOppFaceIndexOfHex
1701 //purpose  : Return index of the opposite face of the hexahedron
1702 //=======================================================================
1703
1704 int SMDS_VolumeTool::GetOppFaceIndexOfHex( int faceIndex )
1705 {
1706   return Hexa_oppF[ faceIndex ];
1707 }
1708
1709 //=======================================================================
1710 //function : IsLinked
1711 //purpose  : return true if theNode1 is linked with theNode2
1712 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1713 //=======================================================================
1714
1715 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
1716                                 const SMDS_MeshNode* theNode2,
1717                                 const bool           theIgnoreMediumNodes) const
1718 {
1719   if ( !myVolume )
1720     return false;
1721
1722   if (myVolume->IsPoly()) {
1723     if (!myPolyedre) {
1724       MESSAGE("Warning: bad volumic element");
1725       return false;
1726     }
1727     if ( !myAllFacesNbNodes ) {
1728       SMDS_VolumeTool*  me = const_cast< SMDS_VolumeTool* >( this );
1729       me->myPolyQuantities = myPolyedre->GetQuantities();
1730       myAllFacesNbNodes    = &myPolyQuantities[0];
1731     }
1732     int from, to = 0, d1 = 1, d2 = 2;
1733     if ( myPolyedre->IsQuadratic() ) {
1734       if ( theIgnoreMediumNodes ) {
1735         d1 = 2; d2 = 0;
1736       }
1737     } else {
1738       d2 = 0;
1739     }
1740     std::vector<const SMDS_MeshNode*>::const_iterator i;
1741     for (int iface = 0; iface < myNbFaces; iface++)
1742     {
1743       from = to;
1744       to  += myPolyQuantities[iface];
1745       i    = std::find( myVolumeNodes.begin() + from, myVolumeNodes.begin() + to, theNode1 );
1746       if ( i != myVolumeNodes.end() )
1747       {
1748         if ((  theNode2 == *( i-d1 ) ||
1749                theNode2 == *( i+d1 )))
1750           return true;
1751         if (( d2 ) &&
1752             (( theNode2 == *( i-d2 ) ||
1753                theNode2 == *( i+d2 ))))
1754           return true;
1755       }
1756     }
1757     return false;
1758   }
1759
1760   // find nodes indices
1761   int i1 = -1, i2 = -1, nbFound = 0;
1762   for ( size_t i = 0; i < myVolumeNodes.size() && nbFound < 2; i++ )
1763   {
1764     if ( myVolumeNodes[ i ] == theNode1 )
1765       i1 = i, ++nbFound;
1766     else if ( myVolumeNodes[ i ] == theNode2 )
1767       i2 = i, ++nbFound;
1768   }
1769   return IsLinked( i1, i2 );
1770 }
1771
1772 //=======================================================================
1773 //function : IsLinked
1774 //purpose  : return true if the node with theNode1Index is linked
1775 //           with the node with theNode2Index
1776 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1777 //=======================================================================
1778
1779 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
1780                                 const int theNode2Index,
1781                                 bool      theIgnoreMediumNodes) const
1782 {
1783   if ( myVolume->IsPoly() ) {
1784     return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
1785   }
1786
1787   int minInd = std::min( theNode1Index, theNode2Index );
1788   int maxInd = std::max( theNode1Index, theNode2Index );
1789
1790   if ( minInd < 0 || maxInd > (int)myVolumeNodes.size() - 1 || maxInd == minInd )
1791     return false;
1792
1793   VolumeType type = GetVolumeType();
1794   if ( myVolume->IsQuadratic() )
1795   {
1796     int firstMediumInd = myVolume->NbCornerNodes();
1797     if ( minInd >= firstMediumInd )
1798       return false; // both nodes are medium - not linked
1799     if ( maxInd < firstMediumInd ) // both nodes are corners
1800     {
1801       if ( theIgnoreMediumNodes )
1802         type = quadToLinear(type); // to check linkage of corner nodes only
1803       else
1804         return false; // corner nodes are not linked directly in a quadratic cell
1805     }
1806   }
1807
1808   switch ( type ) {
1809   case TETRA:
1810     return true;
1811   case HEXA:
1812     switch ( maxInd - minInd ) {
1813     case 1: return minInd != 3;
1814     case 3: return minInd == 0 || minInd == 4;
1815     case 4: return true;
1816     default:;
1817     }
1818     break;
1819   case PYRAM:
1820     if ( maxInd == 4 )
1821       return true;
1822     switch ( maxInd - minInd ) {
1823     case 1:
1824     case 3: return true;
1825     default:;
1826     }
1827     break;
1828   case PENTA:
1829     switch ( maxInd - minInd ) {
1830     case 1: return minInd != 2;
1831     case 2: return minInd == 0 || minInd == 3;
1832     case 3: return true;
1833     default:;
1834     }
1835     break;
1836   case QUAD_TETRA:
1837     {
1838       switch ( minInd ) {
1839       case 0: return ( maxInd==4 ||  maxInd==6 ||  maxInd==7 );
1840       case 1: return ( maxInd==4 ||  maxInd==5 ||  maxInd==8 );
1841       case 2: return ( maxInd==5 ||  maxInd==6 ||  maxInd==9 );
1842       case 3: return ( maxInd==7 ||  maxInd==8 ||  maxInd==9 );
1843       default:;
1844       }
1845       break;
1846     }
1847   case QUAD_HEXA:
1848     {
1849       switch ( minInd ) {
1850       case 0: return ( maxInd==8  ||  maxInd==11 ||  maxInd==16 );
1851       case 1: return ( maxInd==8  ||  maxInd==9  ||  maxInd==17 );
1852       case 2: return ( maxInd==9  ||  maxInd==10 ||  maxInd==18 );
1853       case 3: return ( maxInd==10 ||  maxInd==11 ||  maxInd==19 );
1854       case 4: return ( maxInd==12 ||  maxInd==15 ||  maxInd==16 );
1855       case 5: return ( maxInd==12 ||  maxInd==13 ||  maxInd==17 );
1856       case 6: return ( maxInd==13 ||  maxInd==14 ||  maxInd==18 );
1857       case 7: return ( maxInd==14 ||  maxInd==15 ||  maxInd==19 );
1858       default:;
1859       }
1860       break;
1861     }
1862   case QUAD_PYRAM:
1863     {
1864       switch ( minInd ) {
1865       case 0: return ( maxInd==5 ||  maxInd==8  ||  maxInd==9  );
1866       case 1: return ( maxInd==5 ||  maxInd==6  ||  maxInd==10 );
1867       case 2: return ( maxInd==6 ||  maxInd==7  ||  maxInd==11 );
1868       case 3: return ( maxInd==7 ||  maxInd==8  ||  maxInd==12 );
1869       case 4: return ( maxInd==9 ||  maxInd==10 ||  maxInd==11 ||  maxInd==12 );
1870       default:;
1871       }
1872       break;
1873     }
1874   case QUAD_PENTA:
1875     {
1876       switch ( minInd ) {
1877       case 0: return ( maxInd==6  ||  maxInd==8  ||  maxInd==12 );
1878       case 1: return ( maxInd==6  ||  maxInd==7  ||  maxInd==13 );
1879       case 2: return ( maxInd==7  ||  maxInd==8  ||  maxInd==14 );
1880       case 3: return ( maxInd==9  ||  maxInd==11 ||  maxInd==12 );
1881       case 4: return ( maxInd==9  ||  maxInd==10 ||  maxInd==13 );
1882       case 5: return ( maxInd==10 ||  maxInd==11 ||  maxInd==14 );
1883       default:;
1884       }
1885       break;
1886     }
1887   case HEX_PRISM:
1888     {
1889       const int diff = maxInd-minInd;
1890       if ( diff > 6  ) return false;// not linked top and bottom
1891       if ( diff == 6 ) return true; // linked top and bottom
1892       return diff == 1 || diff == 7;
1893     }
1894   default:;
1895   }
1896   return false;
1897 }
1898
1899 //=======================================================================
1900 //function : GetNodeIndex
1901 //purpose  : Return an index of theNode
1902 //=======================================================================
1903
1904 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
1905 {
1906   if ( myVolume ) {
1907     for ( size_t i = 0; i < myVolumeNodes.size(); i++ ) {
1908       if ( myVolumeNodes[ i ] == theNode )
1909         return i;
1910     }
1911   }
1912   return -1;
1913 }
1914
1915 //================================================================================
1916 /*!
1917  * \brief Fill vector with boundary faces existing in the mesh
1918   * \param faces - vector of found nodes
1919   * \retval int - nb of found faces
1920  */
1921 //================================================================================
1922
1923 int SMDS_VolumeTool::GetAllExistingFaces(std::vector<const SMDS_MeshElement*> & faces) const
1924 {
1925   faces.clear();
1926   SaveFacet savedFacet( myCurFace );
1927   if ( IsPoly() )
1928     for ( int iF = 0; iF < NbFaces(); ++iF ) {
1929       if ( setFace( iF ))
1930         if ( const SMDS_MeshElement* face = SMDS_Mesh::FindFace( myCurFace.myNodes ))
1931           faces.push_back( face );
1932     }
1933   else
1934     for ( int iF = 0; iF < NbFaces(); ++iF ) {
1935       const SMDS_MeshFace* face = 0;
1936       const SMDS_MeshNode** nodes = GetFaceNodes( iF );
1937       switch ( NbFaceNodes( iF )) {
1938       case 3:
1939         face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2] ); break;
1940       case 4:
1941         face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
1942       case 6:
1943         face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2],
1944                                     nodes[3], nodes[4], nodes[5]); break;
1945       case 8:
1946         face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3],
1947                                     nodes[4], nodes[5], nodes[6], nodes[7]); break;
1948       }
1949       if ( face )
1950         faces.push_back( face );
1951     }
1952   return faces.size();
1953 }
1954
1955
1956 //================================================================================
1957 /*!
1958  * \brief Fill vector with boundary edges existing in the mesh
1959  * \param edges - vector of found edges
1960  * \retval int - nb of found faces
1961  */
1962 //================================================================================
1963
1964 int SMDS_VolumeTool::GetAllExistingEdges(std::vector<const SMDS_MeshElement*> & edges) const
1965 {
1966   edges.clear();
1967   edges.reserve( myVolumeNodes.size() * 2 );
1968   for ( size_t i = 0; i < myVolumeNodes.size()-1; ++i ) {
1969     for ( size_t j = i + 1; j < myVolumeNodes.size(); ++j ) {
1970       if ( IsLinked( i, j )) {
1971         const SMDS_MeshElement* edge =
1972           SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] );
1973         if ( edge )
1974           edges.push_back( edge );
1975       }
1976     }
1977   }
1978   return edges.size();
1979 }
1980
1981 //================================================================================
1982 /*!
1983  * \brief Return minimal square distance between connected corner nodes
1984  */
1985 //================================================================================
1986
1987 double SMDS_VolumeTool::MinLinearSize2() const
1988 {
1989   double minSize = 1e+100;
1990   int iQ = myVolume->IsQuadratic() ? 2 : 1;
1991
1992   SaveFacet savedFacet( myCurFace );
1993
1994   // it seems that compute distance twice is faster than organization of a sole computing
1995   myCurFace.myIndex = -1;
1996   for ( int iF = 0; iF < myNbFaces; ++iF )
1997   {
1998     setFace( iF );
1999     for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
2000     {
2001       XYZ n1( myCurFace.myNodes[ iN ]);
2002       XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
2003       minSize = std::min( minSize, (n1 - n2).SquareMagnitude());
2004     }
2005   }
2006
2007   return minSize;
2008 }
2009
2010 //================================================================================
2011 /*!
2012  * \brief Return maximal square distance between connected corner nodes
2013  */
2014 //================================================================================
2015
2016 double SMDS_VolumeTool::MaxLinearSize2() const
2017 {
2018   double maxSize = -1e+100;
2019   int iQ = myVolume->IsQuadratic() ? 2 : 1;
2020
2021   SaveFacet savedFacet( myCurFace );
2022
2023   // it seems that compute distance twice is faster than organization of a sole computing
2024   myCurFace.myIndex = -1;
2025   for ( int iF = 0; iF < myNbFaces; ++iF )
2026   {
2027     setFace( iF );
2028     for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
2029     {
2030       XYZ n1( myCurFace.myNodes[ iN ]);
2031       XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
2032       maxSize = std::max( maxSize, (n1 - n2).SquareMagnitude());
2033     }
2034   }
2035
2036   return maxSize;
2037 }
2038
2039 //================================================================================
2040 /*!
2041  * \brief Fast quickly check that only one volume is built on the face nodes
2042  *        This check is valid for conformal meshes only
2043  */
2044 //================================================================================
2045
2046 bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
2047 {
2048   const bool isFree = true;
2049
2050   if ( !setFace( faceIndex ))
2051     return !isFree;
2052
2053   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
2054
2055   const int  di = myVolume->IsQuadratic() ? 2 : 1;
2056   const int nbN = ( myCurFace.myNbNodes/di <= 4 && !IsPoly()) ? 3 : myCurFace.myNbNodes/di; // nb nodes to check
2057
2058   SMDS_ElemIteratorPtr eIt = nodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
2059   while ( eIt->more() )
2060   {
2061     const SMDS_MeshElement* vol = eIt->next();
2062     if ( vol == myVolume )
2063       continue;
2064     int iN;
2065     for ( iN = 1; iN < nbN; ++iN )
2066       if ( vol->GetNodeIndex( nodes[ iN*di ]) < 0 )
2067         break;
2068     if ( iN == nbN ) // nbN nodes are shared with vol
2069     {
2070       // if ( vol->IsPoly() || vol->NbFaces() > 6 ) // vol is polyhed or hex prism
2071       // {
2072       //   int nb = myCurFace.myNbNodes;
2073       //   if ( myVolume->GetEntityType() != vol->GetEntityType() )
2074       //     nb -= ( GetCenterNodeIndex(0) > 0 );
2075       //   std::set<const SMDS_MeshNode*> faceNodes( nodes, nodes + nb );
2076       //   if ( SMDS_VolumeTool( vol ).GetFaceIndex( faceNodes ) < 0 )
2077       //     continue;
2078       // }
2079       if ( otherVol ) *otherVol = vol;
2080       return !isFree;
2081     }
2082   }
2083   if ( otherVol ) *otherVol = 0;
2084   return isFree;
2085 }
2086
2087 //================================================================================
2088 /*!
2089  * \brief Thorough check that only one volume is built on the face nodes
2090  */
2091 //================================================================================
2092
2093 bool SMDS_VolumeTool::IsFreeFaceAdv( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
2094 {
2095   const bool isFree = true;
2096
2097   if (!setFace( faceIndex ))
2098     return !isFree;
2099
2100   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
2101   const int nbFaceNodes = myCurFace.myNbNodes;
2102
2103   // evaluate nb of face nodes shared by other volumes
2104   int maxNbShared = -1;
2105   typedef std::map< const SMDS_MeshElement*, int > TElemIntMap;
2106   TElemIntMap volNbShared;
2107   TElemIntMap::iterator vNbIt;
2108   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
2109     const SMDS_MeshNode* n = nodes[ iNode ];
2110     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
2111     while ( eIt->more() ) {
2112       const SMDS_MeshElement* elem = eIt->next();
2113       if ( elem != myVolume ) {
2114         vNbIt = volNbShared.insert( std::make_pair( elem, 0 )).first;
2115         (*vNbIt).second++;
2116         if ( vNbIt->second > maxNbShared )
2117           maxNbShared = vNbIt->second;
2118       }
2119     }
2120   }
2121   if ( maxNbShared < 3 )
2122     return isFree; // is free
2123
2124   // find volumes laying on the opposite side of the face
2125   // and sharing all nodes
2126   XYZ intNormal; // internal normal
2127   GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
2128   if ( IsFaceExternal( faceIndex ))
2129     intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
2130   XYZ p0 ( nodes[0] ), baryCenter;
2131   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end();  ) {
2132     const int& nbShared = (*vNbIt).second;
2133     if ( nbShared >= 3 ) {
2134       SMDS_VolumeTool volume( (*vNbIt).first );
2135       volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
2136       XYZ intNormal2( baryCenter - p0 );
2137       if ( intNormal.Dot( intNormal2 ) < 0 ) {
2138         // opposite side
2139         if ( nbShared >= nbFaceNodes )
2140         {
2141           // a volume shares the whole facet
2142           if ( otherVol ) *otherVol = vNbIt->first;
2143           return !isFree;
2144         }
2145         ++vNbIt;
2146         continue;
2147       }
2148     }
2149     // remove a volume from volNbShared map
2150     volNbShared.erase( vNbIt++ );
2151   }
2152
2153   // here volNbShared contains only volumes laying on the opposite side of
2154   // the face and sharing 3 or more but not all face nodes with myVolume
2155   if ( volNbShared.size() < 2 ) {
2156     return isFree; // is free
2157   }
2158
2159   // check if the whole area of a face is shared
2160   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
2161   {
2162     const SMDS_MeshNode* n = nodes[ iNode ];
2163     // check if n is shared by one of volumes of volNbShared
2164     bool isShared = false;
2165     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
2166     while ( eIt->more() && !isShared )
2167       isShared = volNbShared.count( eIt->next() );
2168     if ( !isShared )
2169       return isFree;
2170   }
2171   if ( otherVol ) *otherVol = volNbShared.begin()->first;
2172   return !isFree;
2173
2174 //   if ( !myVolume->IsPoly() )
2175 //   {
2176 //     bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
2177 //     for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
2178 //       SMDS_VolumeTool volume( (*vNbIt).first );
2179 //       bool prevLinkShared = false;
2180 //       int nbSharedLinks = 0;
2181 //       for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
2182 //         bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
2183 //         if ( linkShared )
2184 //           nbSharedLinks++;
2185 //         if ( linkShared && prevLinkShared &&
2186 //              volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
2187 //           isShared[ iNode ] = true;
2188 //         prevLinkShared = linkShared;
2189 //       }
2190 //       if ( nbSharedLinks == nbFaceNodes )
2191 //         return !free; // is not free
2192 //       if ( nbFaceNodes == 4 ) {
2193 //         // check traingle parts 1 & 3
2194 //         if ( isShared[1] && isShared[3] )
2195 //           return !free; // is not free
2196 //         // check triangle parts 0 & 2;
2197 //         // 0 part could not be checked in the loop; check it here
2198 //         if ( isShared[2] && prevLinkShared &&
2199 //              volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
2200 //              volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
2201 //           return !free; // is not free
2202 //       }
2203 //     }
2204 //   }
2205 //  return free;
2206 }
2207
2208 //=======================================================================
2209 //function : GetFaceIndex
2210 //purpose  : Return index of a face formed by theFaceNodes
2211 //=======================================================================
2212
2213 int SMDS_VolumeTool::GetFaceIndex( const std::set<const SMDS_MeshNode*>& theFaceNodes,
2214                                    const int                        theFaceIndexHint ) const
2215 {
2216   if ( theFaceIndexHint >= 0 )
2217   {
2218     int nbNodes = NbFaceNodes( theFaceIndexHint );
2219     if ( nbNodes == (int) theFaceNodes.size() )
2220     {
2221       const SMDS_MeshNode** nodes = GetFaceNodes( theFaceIndexHint );
2222       while ( nbNodes )
2223         if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
2224           --nbNodes;
2225         else
2226           break;
2227       if ( nbNodes == 0 )
2228         return theFaceIndexHint;
2229     }
2230   }
2231   for ( int iFace = 0; iFace < myNbFaces; iFace++ )
2232   {
2233     if ( iFace == theFaceIndexHint )
2234       continue;
2235     int nbNodes = NbFaceNodes( iFace );
2236     if ( nbNodes == (int) theFaceNodes.size() )
2237     {
2238       const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
2239       while ( nbNodes )
2240         if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
2241           --nbNodes;
2242         else
2243           break;
2244       if ( nbNodes == 0 )
2245         return iFace;
2246     }
2247   }
2248   return -1;
2249 }
2250
2251 //=======================================================================
2252 //function : GetFaceIndex
2253 //purpose  : Return index of a face formed by theFaceNodes
2254 //=======================================================================
2255
2256 /*int SMDS_VolumeTool::GetFaceIndex( const std::set<int>& theFaceNodesIndices )
2257 {
2258   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
2259     const int* nodes = GetFaceNodesIndices( iFace );
2260     int nbFaceNodes = NbFaceNodes( iFace );
2261     std::set<int> nodeSet;
2262     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
2263       nodeSet.insert( nodes[ iNode ] );
2264     if ( theFaceNodesIndices == nodeSet )
2265       return iFace;
2266   }
2267   return -1;
2268 }*/
2269
2270 //=======================================================================
2271 //function : setFace
2272 //purpose  :
2273 //=======================================================================
2274
2275 bool SMDS_VolumeTool::setFace( int faceIndex ) const
2276 {
2277   if ( !myVolume )
2278     return false;
2279
2280   if ( myCurFace.myIndex == faceIndex )
2281     return true;
2282
2283   myCurFace.myIndex = -1;
2284
2285   if ( faceIndex < 0 || faceIndex >= NbFaces() )
2286     return false;
2287
2288   if (myVolume->IsPoly())
2289   {
2290     if ( !myPolyedre ) {
2291       MESSAGE("Warning: bad volumic element");
2292       return false;
2293     }
2294
2295     // set face nodes
2296     SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
2297     if ( !myAllFacesNbNodes ) {
2298       me->myPolyQuantities = myPolyedre->GetQuantities();
2299       myAllFacesNbNodes    = &myPolyQuantities[0];
2300     }
2301     myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
2302     myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
2303     me->myPolyIndices.resize( myCurFace.myNbNodes + 1 );
2304     myCurFace.myNodeIndices = & me->myPolyIndices[0];
2305     int shift = std::accumulate( myAllFacesNbNodes, myAllFacesNbNodes+faceIndex, 0 );
2306     for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
2307     {
2308       myCurFace.myNodes      [ iNode ] = myVolumeNodes[ shift + iNode ];
2309       myCurFace.myNodeIndices[ iNode ] = shift + iNode;
2310     }
2311     myCurFace.myNodes      [ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ]; // last = first
2312     myCurFace.myNodeIndices[ myCurFace.myNbNodes ] = myCurFace.myNodeIndices[ 0 ];
2313
2314     // check orientation
2315     if (myExternalFaces)
2316     {
2317       myCurFace.myIndex = faceIndex; // avoid infinite recursion in IsFaceExternal()
2318       myExternalFaces = false; // force normal computation by IsFaceExternal()
2319       if ( !IsFaceExternal( faceIndex ))
2320         std::reverse( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
2321       myExternalFaces = true;
2322     }
2323   }
2324   else
2325   {
2326     if ( !myAllFacesNodeIndices_F )
2327     {
2328       // choose data for an element type
2329       switch ( myVolumeNodes.size() ) {
2330       case 4:
2331         myAllFacesNodeIndices_F  = &Tetra_F [0][0];
2332         //myAllFacesNodeIndices_FE = &Tetra_F [0][0];
2333         myAllFacesNodeIndices_RE = &Tetra_RE[0][0];
2334         myAllFacesNbNodes        = Tetra_nbN;
2335         myMaxFaceNbNodes         = sizeof(Tetra_F[0])/sizeof(Tetra_F[0][0]);
2336         break;
2337       case 5:
2338         myAllFacesNodeIndices_F  = &Pyramid_F [0][0];
2339         //myAllFacesNodeIndices_FE = &Pyramid_F [0][0];
2340         myAllFacesNodeIndices_RE = &Pyramid_RE[0][0];
2341         myAllFacesNbNodes        = Pyramid_nbN;
2342         myMaxFaceNbNodes         = sizeof(Pyramid_F[0])/sizeof(Pyramid_F[0][0]);
2343         break;
2344       case 6:
2345         myAllFacesNodeIndices_F  = &Penta_F [0][0];
2346         //myAllFacesNodeIndices_FE = &Penta_FE[0][0];
2347         myAllFacesNodeIndices_RE = &Penta_RE[0][0];
2348         myAllFacesNbNodes        = Penta_nbN;
2349         myMaxFaceNbNodes         = sizeof(Penta_F[0])/sizeof(Penta_F[0][0]);
2350         break;
2351       case 8:
2352         myAllFacesNodeIndices_F  = &Hexa_F [0][0];
2353         ///myAllFacesNodeIndices_FE = &Hexa_FE[0][0];
2354         myAllFacesNodeIndices_RE = &Hexa_RE[0][0];
2355         myAllFacesNbNodes        = Hexa_nbN;
2356         myMaxFaceNbNodes         = sizeof(Hexa_F[0])/sizeof(Hexa_F[0][0]);
2357         break;
2358       case 10:
2359         myAllFacesNodeIndices_F  = &QuadTetra_F [0][0];
2360         //myAllFacesNodeIndices_FE = &QuadTetra_F [0][0];
2361         myAllFacesNodeIndices_RE = &QuadTetra_RE[0][0];
2362         myAllFacesNbNodes        = QuadTetra_nbN;
2363         myMaxFaceNbNodes         = sizeof(QuadTetra_F[0])/sizeof(QuadTetra_F[0][0]);
2364         break;
2365       case 13:
2366         myAllFacesNodeIndices_F  = &QuadPyram_F [0][0];
2367         //myAllFacesNodeIndices_FE = &QuadPyram_F [0][0];
2368         myAllFacesNodeIndices_RE = &QuadPyram_RE[0][0];
2369         myAllFacesNbNodes        = QuadPyram_nbN;
2370         myMaxFaceNbNodes         = sizeof(QuadPyram_F[0])/sizeof(QuadPyram_F[0][0]);
2371         break;
2372       case 15:
2373         myAllFacesNodeIndices_F  = &QuadPenta_F [0][0];
2374         //myAllFacesNodeIndices_FE = &QuadPenta_FE[0][0];
2375         myAllFacesNodeIndices_RE = &QuadPenta_RE[0][0];
2376         myAllFacesNbNodes        = QuadPenta_nbN;
2377         myMaxFaceNbNodes         = sizeof(QuadPenta_F[0])/sizeof(QuadPenta_F[0][0]);
2378         break;
2379       case 20:
2380       case 27:
2381         myAllFacesNodeIndices_F  = &QuadHexa_F [0][0];
2382         //myAllFacesNodeIndices_FE = &QuadHexa_FE[0][0];
2383         myAllFacesNodeIndices_RE = &QuadHexa_RE[0][0];
2384         myAllFacesNbNodes        = QuadHexa_nbN;
2385         myMaxFaceNbNodes         = sizeof(QuadHexa_F[0])/sizeof(QuadHexa_F[0][0]);
2386         if ( !myIgnoreCentralNodes && myVolumeNodes.size() == 27 )
2387         {
2388           myAllFacesNodeIndices_F  = &TriQuadHexa_F [0][0];
2389           //myAllFacesNodeIndices_FE = &TriQuadHexa_FE[0][0];
2390           myAllFacesNodeIndices_RE = &TriQuadHexa_RE[0][0];
2391           myAllFacesNbNodes        = TriQuadHexa_nbN;
2392           myMaxFaceNbNodes         = sizeof(TriQuadHexa_F[0])/sizeof(TriQuadHexa_F[0][0]);
2393         }
2394         break;
2395       case 12:
2396         myAllFacesNodeIndices_F  = &HexPrism_F [0][0];
2397         //myAllFacesNodeIndices_FE = &HexPrism_FE[0][0];
2398         myAllFacesNodeIndices_RE = &HexPrism_RE[0][0];
2399         myAllFacesNbNodes        = HexPrism_nbN;
2400         myMaxFaceNbNodes         = sizeof(HexPrism_F[0])/sizeof(HexPrism_F[0][0]);
2401         break;
2402       default:
2403         return false;
2404       }
2405     }
2406     myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
2407     // if ( myExternalFaces )
2408     //   myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2409     // else
2410     //   myCurFace.myNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes );
2411     myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
2412
2413     // set face nodes
2414     myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
2415     for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
2416       myCurFace.myNodes[ iNode ] = myVolumeNodes[ myCurFace.myNodeIndices[ iNode ]];
2417     myCurFace.myNodes[ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ];
2418   }
2419
2420   myCurFace.myIndex = faceIndex;
2421
2422   return true;
2423 }
2424
2425 //=======================================================================
2426 //function : GetType
2427 //purpose  : return VolumeType by nb of nodes in a volume
2428 //=======================================================================
2429
2430 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
2431 {
2432   switch ( nbNodes ) {
2433   case 4: return TETRA;
2434   case 5: return PYRAM;
2435   case 6: return PENTA;
2436   case 8: return HEXA;
2437   case 10: return QUAD_TETRA;
2438   case 13: return QUAD_PYRAM;
2439   case 15: return QUAD_PENTA;
2440   case 20:
2441   case 27: return QUAD_HEXA;
2442   case 12: return HEX_PRISM;
2443   default:return UNKNOWN;
2444   }
2445 }
2446
2447 //=======================================================================
2448 //function : NbFaces
2449 //purpose  : return nb of faces by volume type
2450 //=======================================================================
2451
2452 int SMDS_VolumeTool::NbFaces( VolumeType type )
2453 {
2454   switch ( type ) {
2455   case TETRA     :
2456   case QUAD_TETRA: return 4;
2457   case PYRAM     :
2458   case QUAD_PYRAM: return 5;
2459   case PENTA     :
2460   case QUAD_PENTA: return 5;
2461   case HEXA      :
2462   case QUAD_HEXA : return 6;
2463   case HEX_PRISM : return 8;
2464   default:         return 0;
2465   }
2466 }
2467
2468 //================================================================================
2469 /*!
2470  * \brief Useful to know nb of corner nodes of a quadratic volume
2471   * \param type - volume type
2472   * \retval int - nb of corner nodes
2473  */
2474 //================================================================================
2475
2476 int SMDS_VolumeTool::NbCornerNodes(VolumeType type)
2477 {
2478   switch ( type ) {
2479   case TETRA     :
2480   case QUAD_TETRA: return 4;
2481   case PYRAM     :
2482   case QUAD_PYRAM: return 5;
2483   case PENTA     :
2484   case QUAD_PENTA: return 6;
2485   case HEXA      :
2486   case QUAD_HEXA : return 8;
2487   case HEX_PRISM : return 12;
2488   default:         return 0;
2489   }
2490   return 0;
2491 }
2492   //
2493
2494 //=======================================================================
2495 //function : GetFaceNodesIndices
2496 //purpose  : Return the array of face nodes indices
2497 //           To comfort link iteration, the array
2498 //           length == NbFaceNodes( faceIndex ) + 1 and
2499 //           the last node index == the first one.
2500 //=======================================================================
2501
2502 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
2503                                                 int        faceIndex,
2504                                                 bool       external)
2505 {
2506   switch ( type ) {
2507   case TETRA: return Tetra_F[ faceIndex ];
2508   case PYRAM: return Pyramid_F[ faceIndex ];
2509   case PENTA: return external ? Penta_F[ faceIndex ] : Penta_F[ faceIndex ];
2510   case HEXA:  return external ? Hexa_F[ faceIndex ] : Hexa_F[ faceIndex ];
2511   case QUAD_TETRA: return QuadTetra_F[ faceIndex ];
2512   case QUAD_PYRAM: return QuadPyram_F[ faceIndex ];
2513   case QUAD_PENTA: return external ? QuadPenta_F[ faceIndex ] : QuadPenta_F[ faceIndex ];
2514     // what about SMDSEntity_TriQuad_Hexa?
2515   case QUAD_HEXA:  return external ? QuadHexa_F[ faceIndex ] : QuadHexa_F[ faceIndex ];
2516   case HEX_PRISM:  return external ? HexPrism_F[ faceIndex ] : HexPrism_F[ faceIndex ];
2517   default:;
2518   }
2519   return 0;
2520 }
2521
2522 //=======================================================================
2523 //function : NbFaceNodes
2524 //purpose  : Return number of nodes in the array of face nodes
2525 //=======================================================================
2526
2527 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
2528                                  int        faceIndex )
2529 {
2530   switch ( type ) {
2531   case TETRA: return Tetra_nbN[ faceIndex ];
2532   case PYRAM: return Pyramid_nbN[ faceIndex ];
2533   case PENTA: return Penta_nbN[ faceIndex ];
2534   case HEXA:  return Hexa_nbN[ faceIndex ];
2535   case QUAD_TETRA: return QuadTetra_nbN[ faceIndex ];
2536   case QUAD_PYRAM: return QuadPyram_nbN[ faceIndex ];
2537   case QUAD_PENTA: return QuadPenta_nbN[ faceIndex ];
2538     // what about SMDSEntity_TriQuad_Hexa?
2539   case QUAD_HEXA:  return QuadHexa_nbN[ faceIndex ];
2540   case HEX_PRISM:  return HexPrism_nbN[ faceIndex ];
2541   default:;
2542   }
2543   return 0;
2544 }
2545
2546 //=======================================================================
2547 //function : Element
2548 //purpose  : return element
2549 //=======================================================================
2550
2551 const SMDS_MeshVolume* SMDS_VolumeTool::Element() const
2552 {
2553   return static_cast<const SMDS_MeshVolume*>( myVolume );
2554 }
2555
2556 //=======================================================================
2557 //function : ID
2558 //purpose  : return element ID
2559 //=======================================================================
2560
2561 smIdType SMDS_VolumeTool::ID() const
2562 {
2563   return myVolume ? myVolume->GetID() : 0;
2564 }