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