src/SMDS/SMDS_VolumeTool.cxx

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