src/SMDS/SMDS_VolumeTool.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 //  Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 // File      : SMDS_VolumeTool.cxx
  23 // Created   : Tue Jul 13 12:22:13 2004
  24 // Author    : Edward AGAPOV (eap)
  25 //
  26 #ifdef _MSC_VER
  27 #pragma warning(disable:4786)
  28 #endif
  29
  30 #include "SMDS_VolumeTool.hxx"
  31
  32 #include "SMDS_MeshElement.hxx"
  33 #include "SMDS_MeshNode.hxx"
  34 #include "SMDS_PolyhedralVolumeOfNodes.hxx"
  35 #include "SMDS_Mesh.hxx"
  36
  37 #include "utilities.h"
  38
  39 #include <map>
  40 #include <float.h>
  41 #include <math.h>
  42
  43 using namespace std;
  44
  45 // ======================================================
  46 // Node indices in faces depending on volume orientation
  47 // making most faces normals external
  48 // ======================================================
  49
  50 /*
  51 //           N3
  52 //           +
  53 //          /|\
  54 //         / | \
  55 //        /  |  \
  56 //    N0 +---|---+ N1                TETRAHEDRON
  57 //       \   |   /
  58 //        \  |  /
  59 //         \ | /
  60 //          \|/
  61 //           +
  62 //           N2
  63 */
  64 static int Tetra_F [4][4] = { // FORWARD == EXTERNAL
  65   { 0, 1, 2, 0 },              // All faces have external normals
  66   { 0, 3, 1, 0 },
  67   { 1, 3, 2, 1 },
  68   { 0, 2, 3, 0 }};
  69 static int Tetra_R [4][4] = { // REVERSED
  70   { 0, 1, 2, 0 },             // All faces but a bottom have external normals
  71   { 0, 1, 3, 0 },
  72   { 1, 2, 3, 1 },
  73   { 0, 3, 2, 0 }};
  74 static int Tetra_RE [4][4] = { // REVERSED -> FORWARD (EXTERNAL)
  75   { 0, 2, 1, 0 },              // All faces have external normals
  76   { 0, 1, 3, 0 },
  77   { 1, 2, 3, 1 },
  78   { 0, 3, 2, 0 }};
  79 static int Tetra_nbN [] = { 3, 3, 3, 3 };
  80
  81 //
  82 //     PYRAMID
  83 //
  84 static int Pyramid_F [5][5] = { // FORWARD == EXTERNAL
  85   { 0, 1, 2, 3, 0 },            // All faces have external normals
  86   { 0, 4, 1, 0, 4 },
  87   { 1, 4, 2, 1, 4 },
  88   { 2, 4, 3, 2, 4 },
  89   { 3, 4, 0, 3, 4 }};
  90 static int Pyramid_R [5][5] = { // REVERSED
  91   { 0, 1, 2, 3, 0 },            // All faces but a bottom have external normals
  92   { 0, 1, 4, 0, 4 },
  93   { 1, 2, 4, 1, 4 },
  94   { 2, 3, 4, 2, 4 },
  95   { 3, 0, 4, 3, 4 }};
  96 static int Pyramid_RE [5][5] = { // REVERSED -> FORWARD (EXTERNAL)
  97   { 0, 3, 2, 1, 0 },             // All faces but a bottom have external normals
  98   { 0, 1, 4, 0, 4 },
  99   { 1, 2, 4, 1, 4 },
 100   { 2, 3, 4, 2, 4 },
 101   { 3, 0, 4, 3, 4 }};
 102 static int Pyramid_nbN [] = { 4, 3, 3, 3, 3 };
 103
 104 /*
 105 //            + N4
 106 //           /|\
 107 //          / | \
 108 //         /  |  \
 109 //        /   |   \
 110 //    N3 +---------+ N5
 111 //       |    |    |
 112 //       |    + N1 |
 113 //       |   / \   |                PENTAHEDRON
 114 //       |  /   \  |
 115 //       | /     \ |
 116 //       |/       \|
 117 //    N0 +---------+ N2
 118 */
 119 static int Penta_F [5][5] = { // FORWARD
 120   { 0, 1, 2, 0, 0 },          // Top face has an internal normal, other - external
 121   { 3, 4, 5, 3, 3 },          // 0 is bottom, 1 is top face
 122   { 0, 2, 5, 3, 0 },
 123   { 1, 4, 5, 2, 1 },
 124   { 0, 3, 4, 1, 0 }};
 125 static int Penta_R [5][5] = { // REVERSED
 126   { 0, 1, 2, 0, 0 },          // Bottom face has an internal normal, other - external
 127   { 3, 4, 5, 3, 3 },          // 0 is bottom, 1 is top face
 128   { 0, 3, 5, 2, 0 },
 129   { 1, 2, 5, 4, 1 },
 130   { 0, 1, 4, 3, 0 }};
 131 static int Penta_FE [5][5] = { // FORWARD -> EXTERNAL
 132   { 0, 1, 2, 0, 0 },
 133   { 3, 5, 4, 3, 3 },
 134   { 0, 2, 5, 3, 0 },
 135   { 1, 4, 5, 2, 1 },
 136   { 0, 3, 4, 1, 0 }};
 137 static int Penta_RE [5][5] = { // REVERSED -> EXTERNAL
 138   { 0, 2, 1, 0, 0 },
 139   { 3, 4, 5, 3, 3 },
 140   { 0, 3, 5, 2, 0 },
 141   { 1, 2, 5, 4, 1 },
 142   { 0, 1, 4, 3, 0 }};
 143 static int Penta_nbN [] = { 3, 3, 4, 4, 4 };
 144
 145 /*
 146 //         N5+----------+N6
 147 //          /|         /|
 148 //         / |        / |
 149 //        /  |       /  |
 150 //     N4+----------+N7 |
 151 //       |   |      |   |           HEXAHEDRON
 152 //       |   |      |   |
 153 //       |   |      |   |
 154 //       | N1+------|---+N2
 155 //       |  /       |  /
 156 //       | /        | /
 157 //       |/         |/
 158 //     N0+----------+N3
 159 */
 160 static int Hexa_F [6][5] = { // FORWARD
 161   { 0, 1, 2, 3, 0 },         // opposite faces are neighbouring,
 162   { 4, 5, 6, 7, 4 },         // odd face(1,3,5) normal is internal, even(0,2,4) - external
 163   { 1, 0, 4, 5, 1 },         // same index nodes of opposite faces are linked
 164   { 2, 3, 7, 6, 2 },
 165   { 0, 3, 7, 4, 0 },
 166   { 1, 2, 6, 5, 1 }};
 167 // static int Hexa_R [6][5] = { // REVERSED
 168 //   { 0, 3, 2, 1, 0 },         // opposite faces are neighbouring,
 169 //   { 4, 7, 6, 5, 4 },         // odd face(1,3,5) normal is external, even(0,2,4) - internal
 170 //   { 1, 5, 4, 0, 1 },         // same index nodes of opposite faces are linked
 171 //   { 2, 6, 7, 3, 2 },
 172 //   { 0, 4, 7, 3, 0 },
 173 //   { 1, 5, 6, 2, 1 }};
 174 static int Hexa_FE [6][5] = { // FORWARD -> EXTERNAL
 175   { 0, 1, 2, 3, 0 } ,         // opposite faces are neighbouring,
 176   { 4, 7, 6, 5, 4 },          // all face normals are external,
 177   { 0, 4, 5, 1, 0 },          // links in opposite faces: 0-0, 1-3, 2-2, 3-1
 178   { 3, 2, 6, 7, 3 },
 179   { 0, 3, 7, 4, 0 },
 180   { 1, 5, 6, 2, 1 }};
 181 static int Hexa_RE [6][5] = { // REVERSED -> EXTERNAL
 182   { 0, 3, 2, 1, 0 },          // opposite faces are neighbouring,
 183   { 4, 5, 6, 7, 4 },          // all face normals are external,
 184   { 0, 1, 5, 4, 0 },          // links in opposite faces: 0-0, 1-3, 2-2, 3-1
 185   { 3, 7, 6, 2, 3 },
 186   { 0, 4, 7, 3, 0 },
 187   { 1, 2, 6, 5, 1 }};
 188 static int Hexa_nbN [] = { 4, 4, 4, 4, 4, 4 };
 189
 190
 191 /*
 192 //           N3
 193 //           +
 194 //          /|\
 195 //        7/ | \8
 196 //        /  |4 \                    QUADRATIC
 197 //    N0 +---|---+ N1                TETRAHEDRON
 198 //       \   +9  /
 199 //        \  |  /
 200 //        6\ | /5
 201 //          \|/
 202 //           +
 203 //           N2
 204 */
 205 static int QuadTetra_F [4][7] = { // FORWARD == EXTERNAL
 206   { 0, 4, 1, 5, 2, 6, 0 },        // All faces have external normals
 207   { 0, 7, 3, 8, 1, 4, 0 },
 208   { 1, 8, 3, 9, 2, 5, 1 },
 209   { 0, 6, 2, 9, 3, 7, 0 }};
 210 static int QuadTetra_R [4][7] = { // REVERSED
 211   { 0, 4, 1, 5, 2, 6, 0 },        // All faces but a bottom have external normals
 212   { 0, 4, 1, 8, 3, 7, 0 },
 213   { 1, 5, 2, 9, 3, 8, 1 },
 214   { 0, 7, 3, 9, 2, 6, 0 }};
 215 static int QuadTetra_RE [4][7] = { // REVERSED -> FORWARD (EXTERNAL)
 216   { 0, 6, 2, 5, 1, 4, 0 },              // All faces have external normals
 217   { 0, 4, 1, 8, 3, 7, 0 },
 218   { 1, 5, 2, 9, 3, 8, 1 },
 219   { 0, 7, 3, 9, 2, 6, 0 }};
 220 static int QuadTetra_nbN [] = { 6, 6, 6, 6 };
 221
 222 //
 223 //     QUADRATIC
 224 //     PYRAMID
 225 //
 226 //            +4
 227 //
 228 //
 229 //       10+-----+11
 230 //         |     |        9 - middle point for (0,4) etc.
 231 //         |     |
 232 //        9+-----+12
 233 //
 234 //            6
 235 //      1+----+----+2
 236 //       |         |
 237 //       |         |
 238 //      5+         +7
 239 //       |         |
 240 //       |         |
 241 //      0+----+----+3
 242 //            8
 243 static int QuadPyram_F [5][9] = {  // FORWARD == EXTERNAL
 244   { 0, 5, 1, 6, 2, 7, 3, 8, 0 },   // All faces have external normals
 245   { 0, 9, 4, 10,1, 5, 0, 4, 4 },
 246   { 1, 10,4, 11,2, 6, 1, 4, 4 },
 247   { 2, 11,4, 12,3, 7, 2, 4, 4 },
 248   { 3, 12,4, 9, 0, 8, 3, 4, 4 }};
 249 static int QuadPyram_R [5][9] = {  // REVERSED
 250   { 0, 5, 1, 6, 2, 7, 3, 8, 0 },   // All faces but a bottom have external normals
 251   { 0, 5, 1, 10,4, 9, 0, 4, 4 },
 252   { 1, 6, 2, 11,4, 10,1, 4, 4 },
 253   { 2, 7, 3, 12,4, 11,2, 4, 4 },
 254   { 3, 8, 0, 9, 4, 12,3, 4, 4 }};
 255 static int QuadPyram_RE [5][9] = { // REVERSED -> FORWARD (EXTERNAL)
 256   { 0, 8, 3, 7, 2, 6, 1, 5, 0 },   // All faces but a bottom have external normals
 257   { 0, 5, 1, 10,4, 9, 0, 4, 4 },
 258   { 1, 6, 2, 11,4, 10,1, 4, 4 },
 259   { 2, 7, 3, 12,4, 11,2, 4, 4 },
 260   { 3, 8, 0, 9, 4, 12,3, 4, 4 }};
 261 static int QuadPyram_nbN [] = { 8, 6, 6, 6, 6 };
 262
 263 /*
 264 //            + N4
 265 //           /|\
 266 //         9/ | \10
 267 //         /  |  \
 268 //        /   |   \
 269 //    N3 +----+----+ N5
 270 //       |    |11  |
 271 //       |    |    |
 272 //       |    +13  |                QUADRATIC
 273 //       |    |    |                PENTAHEDRON
 274 //       |    |    |
 275 //       |    |    |
 276 //       |    |    |
 277 //     12+    |    +14
 278 //       |    |    |
 279 //       |    |    |
 280 //       |    + N1 |
 281 //       |   / \   |
 282 //       | 6/   \7 |
 283 //       | /     \ |
 284 //       |/       \|
 285 //    N0 +---------+ N2
 286 //            8
 287 */
 288 static int QuadPenta_F [5][9] = {  // FORWARD
 289   { 0, 6, 1, 7, 2, 8, 0, 0, 0 },   // Top face has an internal normal, other - external
 290   { 3, 9, 4, 10,5, 11,3, 3, 3 },   // 0 is bottom, 1 is top face
 291   { 0, 8, 2, 14,5, 11,3, 12,0 },
 292   { 1, 13,4, 10,5, 14,2, 7, 1 },
 293   { 0, 12,3, 9, 4, 13,1, 6, 0 }};
 294 static int QuadPenta_R [5][9] = { // REVERSED
 295   { 0, 6, 1, 7, 2, 8, 0, 0, 0 },  // Bottom face has an internal normal, other - external
 296   { 3, 9, 4, 10,5, 11,3, 3, 3 },  // 0 is bottom, 1 is top face
 297   { 0, 12,3, 11,5, 14,2, 8, 0 },
 298   { 1, 7, 2, 14,5, 10,4, 13,1 },
 299   { 0, 6, 1, 13,4, 9, 3, 12,0 }};
 300 static int QuadPenta_FE [5][9] = { // FORWARD -> EXTERNAL
 301   { 0, 6, 1, 7, 2, 8, 0, 0, 0 },
 302   { 3,11, 5, 10,4, 9, 3, 3, 3 },
 303   { 0, 8, 2, 14,5, 11,3, 12,0 },
 304   { 1, 13,4, 10,5, 14,2, 7, 1 },
 305   { 0, 12,3, 9, 4, 13,1, 6, 0 }};
 306 static int QuadPenta_RE [5][9] = { // REVERSED -> EXTERNAL
 307   { 0, 8, 2, 7, 1, 6, 0, 0, 0 },
 308   { 3, 9, 4, 10,5, 11,3, 3, 3 },
 309   { 0, 12,3, 11,5, 14,2, 8, 0 },
 310   { 1, 7, 2, 14,5, 10,4, 13,1 },
 311   { 0, 6, 1, 13,4, 9, 3, 12,0 }};
 312 static int QuadPenta_nbN [] = { 6, 6, 8, 8, 8 };
 313
 314 /*
 315 //                 13
 316 //         N5+-----+-----+N6
 317 //          /|          /|
 318 //       12+ |       14+ |
 319 //        /  |        /  |
 320 //     N4+-----+-----+N7 |           QUADRATIC
 321 //       |   | 15    |   |           HEXAHEDRON
 322 //       |   |       |   |
 323 //       | 17+       |   +18
 324 //       |   |       |   |
 325 //       |   |       |   |
 326 //       |   |       |   |
 327 //     16+   |       +19 |
 328 //       |   |       |   |
 329 //       |   |     9 |   |
 330 //       | N1+-----+-|---+N2
 331 //       |  /        |  /
 332 //       | +8        | +10
 333 //       |/          |/
 334 //     N0+-----+-----+N3
 335 //             11
 336 */
 337 static int QuadHexa_F [6][9] = {  // FORWARD
 338   { 0, 8, 1, 9, 2, 10,3, 11,0 },  // opposite faces are neighbouring,
 339   { 4, 12,5, 13,6, 14,7, 15,4 },  // odd face(1,3,5) normal is internal, even(0,2,4) - external
 340   { 1, 8, 0, 16,4, 12,5, 17,1 },  // same index nodes of opposite faces are linked
 341   { 2, 10,3, 19,7, 14,6, 18,2 },
 342   { 0, 11,3, 19,7, 15,4, 16,0 },
 343   { 1, 9, 2, 18,6, 13,5, 17,1 }};
 344 // static int Hexa_R [6][5] = { // REVERSED
 345 //   { 0, 3, 2, 1, 0 },         // opposite faces are neighbouring,
 346 //   { 4, 7, 6, 5, 4 },         // odd face(1,3,5) normal is external, even(0,2,4) - internal
 347 //   { 1, 5, 4, 0, 1 },         // same index nodes of opposite faces are linked
 348 //   { 2, 6, 7, 3, 2 },
 349 //   { 0, 4, 7, 3, 0 },
 350 //   { 1, 5, 6, 2, 1 }};
 351 static int QuadHexa_FE [6][9] = {  // FORWARD -> EXTERNAL
 352   { 0, 8, 1, 9, 2, 10,3, 11,0 },   // opposite faces are neighbouring,
 353   { 4, 15,7, 14,6, 13,5, 12,4 },   // all face normals are external,
 354   { 0, 16,4, 12,5, 17,1, 8, 0 },   // links in opposite faces: 0-0, 1-3, 2-2, 3-1
 355   { 3, 10,2, 18,6, 14,7, 19,3 },
 356   { 0, 11,3, 19,7, 15,4, 16,0 },
 357   { 1, 17,5, 13,6, 18,2, 9, 1 }};
 358 static int QuadHexa_RE [6][9] = {  // REVERSED -> EXTERNAL
 359   { 0, 11,3, 10,2, 9, 1, 8, 0 },   // opposite faces are neighbouring,
 360   { 4, 12,5, 13,6, 14,7, 15,4 },   // all face normals are external,
 361   { 0, 8, 1, 17,5, 12,4, 16,0 },   // links in opposite faces: 0-0, 1-3, 2-2, 3-1
 362   { 3, 19,7, 14,6, 18,2, 10,3 },
 363   { 0, 16,4, 15,7, 19,3, 11,0 },
 364   { 1, 9, 2, 18,6, 13,5, 17,1 }};
 365 static int QuadHexa_nbN [] = { 8, 8, 8, 8, 8, 8 };
 366
 367
 368 // ========================================================
 369 // to perform some calculations without linkage to CASCADE
 370 // ========================================================
 371 namespace
 372 {
 373 struct XYZ {
 374   double x;
 375   double y;
 376   double z;
 377   XYZ()                               { x = 0; y = 0; z = 0; }
 378   XYZ( double X, double Y, double Z ) { x = X; y = Y; z = Z; }
 379   XYZ( const XYZ& other )             { x = other.x; y = other.y; z = other.z; }
 380   XYZ( const SMDS_MeshNode* n )       { x = n->X(); y = n->Y(); z = n->Z(); }
 381   inline XYZ operator-( const XYZ& other );
 382   inline XYZ Crossed( const XYZ& other );
 383   inline double Dot( const XYZ& other );
 384   inline double Magnitude();
 385 };
 386 inline XYZ XYZ::operator-( const XYZ& Right ) {
 387   return XYZ(x - Right.x, y - Right.y, z - Right.z);
 388 }
 389 inline XYZ XYZ::Crossed( const XYZ& Right ) {
 390   return XYZ (y * Right.z - z * Right.y,
 391               z * Right.x - x * Right.z,
 392               x * Right.y - y * Right.x);
 393 }
 394 inline double XYZ::Dot( const XYZ& Other ) {
 395   return(x * Other.x + y * Other.y + z * Other.z);
 396 }
 397 inline double XYZ::Magnitude() {
 398   return sqrt (x * x + y * y + z * z);
 399 }
 400 }
 401
 402 //=======================================================================
 403 //function : SMDS_VolumeTool
 404 //purpose  :
 405 //=======================================================================
 406
 407 SMDS_VolumeTool::SMDS_VolumeTool ()
 408      : myVolume( 0 ),
 409        myPolyedre( 0 ),
 410        myVolForward( true ),
 411        myNbFaces( 0 ),
 412        myVolumeNbNodes( 0 ),
 413        myVolumeNodes( NULL ),
 414        myExternalFaces( false ),
 415        myFaceNbNodes( 0 ),
 416        myCurFace( -1 ),
 417        myFaceNodeIndices( NULL ),
 418        myFaceNodes( NULL )
 419 {
 420 }
 421
 422 //=======================================================================
 423 //function : SMDS_VolumeTool
 424 //purpose  :
 425 //=======================================================================
 426
 427 SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume)
 428      : myVolume( 0 ),
 429        myPolyedre( 0 ),
 430        myVolForward( true ),
 431        myNbFaces( 0 ),
 432        myVolumeNbNodes( 0 ),
 433        myVolumeNodes( NULL ),
 434        myExternalFaces( false ),
 435        myFaceNbNodes( 0 ),
 436        myCurFace( -1 ),
 437        myFaceNodeIndices( NULL ),
 438        myFaceNodes( NULL )
 439 {
 440   Set( theVolume );
 441 }
 442
 443 //=======================================================================
 444 //function : SMDS_VolumeTool
 445 //purpose  :
 446 //=======================================================================
 447
 448 SMDS_VolumeTool::~SMDS_VolumeTool()
 449 {
 450   if ( myVolumeNodes != NULL ) delete [] myVolumeNodes;
 451   if ( myFaceNodes != NULL   ) delete [] myFaceNodes;
 452
 453   myFaceNodeIndices = NULL;
 454   myVolumeNodes = myFaceNodes = NULL;
 455 }
 456
 457 //=======================================================================
 458 //function : SetVolume
 459 //purpose  : Set volume to iterate on
 460 //=======================================================================
 461
 462 bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume)
 463 {
 464   myVolume = 0;
 465   myPolyedre = 0;
 466
 467   myVolForward = true;
 468   myNbFaces = 0;
 469   myVolumeNbNodes = 0;
 470   if (myVolumeNodes != NULL) {
 471     delete [] myVolumeNodes;
 472     myVolumeNodes = NULL;
 473   }
 474
 475   myExternalFaces = false;
 476   myFaceNbNodes = 0;
 477
 478   myCurFace = -1;
 479   myFaceNodeIndices = NULL;
 480   if (myFaceNodes != NULL) {
 481     delete [] myFaceNodes;
 482     myFaceNodes = NULL;
 483   }
 484
 485   if ( theVolume && theVolume->GetType() == SMDSAbs_Volume )
 486   {
 487     myVolume = theVolume;
 488
 489     myNbFaces = theVolume->NbFaces();
 490     myVolumeNbNodes = theVolume->NbNodes();
 491
 492     // set volume nodes
 493     int iNode = 0;
 494     myVolumeNodes = new const SMDS_MeshNode* [myVolumeNbNodes];
 495     SMDS_ElemIteratorPtr nodeIt = myVolume->nodesIterator();
 496     while ( nodeIt->more() ) {
 497       myVolumeNodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
 498     }
 499
 500     if (myVolume->IsPoly()) {
 501       myPolyedre = static_cast<const SMDS_PolyhedralVolumeOfNodes*>( myVolume );
 502       if (!myPolyedre) {
 503         MESSAGE("Warning: bad volumic element");
 504         return false;
 505       }
 506     }
 507     else {
 508       switch ( myVolumeNbNodes ) {
 509       case 4:
 510       case 5:
 511       case 6:
 512       case 8:
 513       case 10:
 514       case 13:
 515       case 15:
 516       case 20: {
 517         // define volume orientation
 518         XYZ botNormal;
 519         GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z );
 520         const SMDS_MeshNode* topNode = myVolumeNodes[ myVolumeNbNodes - 1 ];
 521         const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ];
 522         XYZ upDir (topNode->X() - botNode->X(),
 523                    topNode->Y() - botNode->Y(),
 524                    topNode->Z() - botNode->Z() );
 525         myVolForward = ( botNormal.Dot( upDir ) < 0 );
 526         break;
 527       }
 528       default:
 529         break;
 530       }
 531     }
 532   }
 533   return ( myVolume != 0 );
 534 }
 535
 536 //=======================================================================
 537 //function : Inverse
 538 //purpose  : Inverse volume
 539 //=======================================================================
 540
 541 #define SWAP_NODES(nodes,i1,i2)           \
 542 {                                         \
 543   const SMDS_MeshNode* tmp = nodes[ i1 ]; \
 544   nodes[ i1 ] = nodes[ i2 ];              \
 545   nodes[ i2 ] = tmp;                      \
 546 }
 547 void SMDS_VolumeTool::Inverse ()
 548 {
 549   if ( !myVolume ) return;
 550
 551   if (myVolume->IsPoly()) {
 552     MESSAGE("Warning: attempt to inverse polyhedral volume");
 553     return;
 554   }
 555
 556   myVolForward = !myVolForward;
 557   myCurFace = -1;
 558
 559   // inverse top and bottom faces
 560   switch ( myVolumeNbNodes ) {
 561   case 4:
 562     SWAP_NODES( myVolumeNodes, 1, 2 );
 563     break;
 564   case 5:
 565     SWAP_NODES( myVolumeNodes, 1, 3 );
 566     break;
 567   case 6:
 568     SWAP_NODES( myVolumeNodes, 1, 2 );
 569     SWAP_NODES( myVolumeNodes, 4, 5 );
 570     break;
 571   case 8:
 572     SWAP_NODES( myVolumeNodes, 1, 3 );
 573     SWAP_NODES( myVolumeNodes, 5, 7 );
 574     break;
 575
 576   case 10:
 577     SWAP_NODES( myVolumeNodes, 1, 2 );
 578     SWAP_NODES( myVolumeNodes, 4, 6 );
 579     SWAP_NODES( myVolumeNodes, 8, 9 );
 580     break;
 581   case 13:
 582     SWAP_NODES( myVolumeNodes, 1, 3 );
 583     SWAP_NODES( myVolumeNodes, 5, 8 );
 584     SWAP_NODES( myVolumeNodes, 6, 7 );
 585     SWAP_NODES( myVolumeNodes, 10, 12 );
 586     break;
 587   case 15:
 588     SWAP_NODES( myVolumeNodes, 1, 2 );
 589     SWAP_NODES( myVolumeNodes, 4, 5 );
 590     SWAP_NODES( myVolumeNodes, 6, 8 );
 591     SWAP_NODES( myVolumeNodes, 9, 11 );
 592     SWAP_NODES( myVolumeNodes, 13, 14 );
 593     break;
 594   case 20:
 595     SWAP_NODES( myVolumeNodes, 1, 3 );
 596     SWAP_NODES( myVolumeNodes, 5, 7 );
 597     SWAP_NODES( myVolumeNodes, 8, 11 );
 598     SWAP_NODES( myVolumeNodes, 9, 10 );
 599     SWAP_NODES( myVolumeNodes, 12, 15 );
 600     SWAP_NODES( myVolumeNodes, 13, 14 );
 601     SWAP_NODES( myVolumeNodes, 17, 19 );
 602     break;
 603   default:;
 604   }
 605 }
 606
 607 //=======================================================================
 608 //function : GetVolumeType
 609 //purpose  :
 610 //=======================================================================
 611
 612 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetVolumeType() const
 613 {
 614   if ( myPolyedre )
 615     return POLYHEDA;
 616
 617   if ( myVolume ) {
 618 //    static const VolumeType types[] = {
 619 //      TETRA,     // myVolumeNbNodes = 4
 620 //      PYRAM,     // myVolumeNbNodes = 5
 621 //      PENTA,     // myVolumeNbNodes = 6
 622 //      UNKNOWN,   // myVolumeNbNodes = 7
 623 //      HEXA       // myVolumeNbNodes = 8
 624 //    };
 625 //    return types[ myVolumeNbNodes - 4 ];
 626     switch(myVolumeNbNodes) {
 627     case 4: return TETRA; break;
 628     case 5: return PYRAM; break;
 629     case 6: return PENTA; break;
 630     case 8: return HEXA; break;
 631     case 10: return QUAD_TETRA; break;
 632     case 13: return QUAD_PYRAM; break;
 633     case 15: return QUAD_PENTA; break;
 634     case 20: return QUAD_HEXA; break;
 635     default: break;
 636     }
 637   }
 638
 639   return UNKNOWN;
 640 }
 641
 642 //=======================================================================
 643 //function : getTetraVolume
 644 //purpose  :
 645 //=======================================================================
 646
 647 static double getTetraVolume(const SMDS_MeshNode* n1,
 648                              const SMDS_MeshNode* n2,
 649                              const SMDS_MeshNode* n3,
 650                              const SMDS_MeshNode* n4)
 651 {
 652   double X1 = n1->X();
 653   double Y1 = n1->Y();
 654   double Z1 = n1->Z();
 655
 656   double X2 = n2->X();
 657   double Y2 = n2->Y();
 658   double Z2 = n2->Z();
 659
 660   double X3 = n3->X();
 661   double Y3 = n3->Y();
 662   double Z3 = n3->Z();
 663
 664   double X4 = n4->X();
 665   double Y4 = n4->Y();
 666   double Z4 = n4->Z();
 667
 668   double Q1 = -(X1-X2)*(Y3*Z4-Y4*Z3);
 669   double Q2 =  (X1-X3)*(Y2*Z4-Y4*Z2);
 670   double R1 = -(X1-X4)*(Y2*Z3-Y3*Z2);
 671   double R2 = -(X2-X3)*(Y1*Z4-Y4*Z1);
 672   double S1 =  (X2-X4)*(Y1*Z3-Y3*Z1);
 673   double S2 = -(X3-X4)*(Y1*Z2-Y2*Z1);
 674
 675   return (Q1+Q2+R1+R2+S1+S2)/6.0;
 676 }
 677
 678 //=======================================================================
 679 //function : GetSize
 680 //purpose  : Return element volume
 681 //=======================================================================
 682
 683 double SMDS_VolumeTool::GetSize() const
 684 {
 685   double V = 0.;
 686   if ( !myVolume )
 687     return 0.;
 688
 689   if ( myVolume->IsPoly() )
 690   {
 691     if ( !myPolyedre )
 692       return 0.;
 693
 694     // split a polyhedron into tetrahedrons
 695
 696     SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
 697     XYZ baryCenter;
 698     me->GetBaryCenter(baryCenter.x, baryCenter.y, baryCenter.z);
 699     SMDS_MeshNode bcNode ( baryCenter.x, baryCenter.y, baryCenter.z );
 700
 701     for ( int f = 0; f < NbFaces(); ++f )
 702     {
 703       bool externalFace = me->IsFaceExternal( f ); // it calls setFace()
 704       for ( int n = 2; n < myFaceNbNodes; ++n )
 705       {
 706         double Vn = getTetraVolume( myFaceNodes[ 0 ],
 707                                     myFaceNodes[ n-1 ],
 708                                     myFaceNodes[ n ],
 709                                     & bcNode );
 710 ///         cout <<"++++   " << Vn << "   nodes " <<myFaceNodes[ 0 ]->GetID() << " " <<myFaceNodes[ n-1 ]->GetID() << " " <<myFaceNodes[ n ]->GetID() << "        < " << V << endl;
 711         V += externalFace ? -Vn : Vn;
 712       }
 713     }
 714   }
 715   else
 716   {
 717     const static int ind[] = {
 718       0, 1, 3, 6, 11, 19, 32, 46, 66};
 719     const static int vtab[][4] = {
 720       // tetrahedron
 721       { 0, 1, 2, 3 },
 722       // pyramid
 723       { 0, 1, 3, 4 },
 724       { 1, 2, 3, 4 },
 725       // pentahedron
 726       { 0, 1, 2, 3 },
 727       { 1, 5, 3, 4 },
 728       { 1, 5, 2, 3 },
 729       // hexahedron
 730       { 1, 4, 3, 0 },
 731       { 4, 1, 6, 5 },
 732       { 1, 3, 6, 2 },
 733       { 4, 6, 3, 7 },
 734       { 1, 4, 6, 3 },
 735
 736       // quadratic tetrahedron
 737       { 0, 4, 6, 7 },
 738       { 1, 5, 4, 8 },
 739       { 2, 6, 5, 9 },
 740       { 7, 8, 9, 3 },
 741       { 4, 6, 7, 9 },
 742       { 4, 5, 6, 9 },
 743       { 4, 7, 8, 9 },
 744       { 4, 5, 9, 8 },
 745
 746       // quadratic pyramid
 747       { 0, 5, 8, 9 },
 748       { 1, 5,10, 6 },
 749       { 2, 6,11, 7 },
 750       { 3, 7,12, 8 },
 751       { 4, 9,11,10 },
 752       { 4, 9,12,11 },
 753       { 10, 5, 9, 8 },
 754       { 10, 8, 9,12 },
 755       { 10, 8,12, 7 },
 756       { 10, 7,12,11 },
 757       { 10, 7,11, 6 },
 758       { 10, 5, 8, 6 },
 759       { 10, 6, 8, 7 },
 760
 761       // quadratic pentahedron
 762       { 12, 0, 8, 6 },
 763       { 12, 8, 7, 6 },
 764       { 12, 8, 2, 7 },
 765       { 12, 6, 7, 1 },
 766       { 12, 1, 7,13 },
 767       { 12, 7, 2,13 },
 768       { 12, 2,14,13 },
 769
 770       { 12, 3, 9,11 },
 771       { 12,11, 9,10 },
 772       { 12,11,10, 5 },
 773       { 12, 9, 4,10 },
 774       { 12,14, 5,10 },
 775       { 12,14,10, 4 },
 776       { 12,14, 4,13 },
 777
 778       // quadratic hexahedron
 779       { 16, 0,11, 8 },
 780       { 16,11, 9, 8 },
 781       { 16, 8, 9, 1 },
 782       { 16,11, 3,10 },
 783       { 16,11,10, 9 },
 784       { 16,10, 2, 9 },
 785       { 16, 3,19, 2 },
 786       { 16, 2,19,18 },
 787       { 16, 2,18,17 },
 788       { 16, 2,17, 1 },
 789
 790       { 16, 4,12,15 },
 791       { 16,12, 5,13 },
 792       { 16,12,13,15 },
 793       { 16,13, 6,14 },
 794       { 16,13,14,15 },
 795       { 16,14, 7,15 },
 796       { 16, 6, 5,17 },
 797       { 16,18, 6,17 },
 798       { 16,18, 7, 6 },
 799       { 16,18,19, 7 },
 800
 801     };
 802
 803     int type = GetVolumeType();
 804     int n1 = ind[type];
 805     int n2 = ind[type+1];
 806
 807     for (int i = n1; i <  n2; i++) {
 808       V -= getTetraVolume( myVolumeNodes[ vtab[i][0] ],
 809                            myVolumeNodes[ vtab[i][1] ],
 810                            myVolumeNodes[ vtab[i][2] ],
 811                            myVolumeNodes[ vtab[i][3] ]);
 812     }
 813   }
 814   return V;
 815 }
 816
 817 //=======================================================================
 818 //function : GetBaryCenter
 819 //purpose  :
 820 //=======================================================================
 821
 822 bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
 823 {
 824   X = Y = Z = 0.;
 825   if ( !myVolume )
 826     return false;
 827
 828   for ( int i = 0; i < myVolumeNbNodes; i++ ) {
 829     X += myVolumeNodes[ i ]->X();
 830     Y += myVolumeNodes[ i ]->Y();
 831     Z += myVolumeNodes[ i ]->Z();
 832   }
 833   X /= myVolumeNbNodes;
 834   Y /= myVolumeNbNodes;
 835   Z /= myVolumeNbNodes;
 836
 837   return true;
 838 }
 839
 840 //================================================================================
 841 /*!
 842  * \brief Classify a point
 843  *  \param tol - thickness of faces
 844  */
 845 //================================================================================
 846
 847 bool SMDS_VolumeTool::IsOut(double X, double Y, double Z, double tol)
 848 {
 849   // LIMITATION: for convex volumes only
 850   XYZ p( X,Y,Z );
 851   for ( int iF = 0; iF < myNbFaces; ++iF )
 852   {
 853     XYZ faceNormal;
 854     if ( !GetFaceNormal( iF, faceNormal.x, faceNormal.y, faceNormal.z ))
 855       continue;
 856     if ( !IsFaceExternal( iF ))
 857       faceNormal = XYZ() - faceNormal; // reverse
 858
 859     XYZ face2p( p - XYZ( myFaceNodes[0] ));
 860     if ( face2p.Dot( faceNormal ) > tol )
 861       return true;
 862   }
 863   return false;
 864 }
 865
 866 //=======================================================================
 867 //function : SetExternalNormal
 868 //purpose  : Node order will be so that faces normals are external
 869 //=======================================================================
 870
 871 void SMDS_VolumeTool::SetExternalNormal ()
 872 {
 873   myExternalFaces = true;
 874   myCurFace = -1;
 875 }
 876
 877 //=======================================================================
 878 //function : NbFaceNodes
 879 //purpose  : Return number of nodes in the array of face nodes
 880 //=======================================================================
 881
 882 int SMDS_VolumeTool::NbFaceNodes( int faceIndex )
 883 {
 884     if ( !setFace( faceIndex ))
 885       return 0;
 886     return myFaceNbNodes;
 887 }
 888
 889 //=======================================================================
 890 //function : GetFaceNodes
 891 //purpose  : Return pointer to the array of face nodes.
 892 //           To comfort link iteration, the array
 893 //           length == NbFaceNodes( faceIndex ) + 1 and
 894 //           the last node == the first one.
 895 //=======================================================================
 896
 897 const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex )
 898 {
 899   if ( !setFace( faceIndex ))
 900     return 0;
 901   return myFaceNodes;
 902 }
 903
 904 //=======================================================================
 905 //function : GetFaceNodesIndices
 906 //purpose  : Return pointer to the array of face nodes indices
 907 //           To comfort link iteration, the array
 908 //           length == NbFaceNodes( faceIndex ) + 1 and
 909 //           the last node index == the first one.
 910 //=======================================================================
 911
 912 const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex )
 913 {
 914   if ( !setFace( faceIndex ))
 915     return 0;
 916
 917   if (myVolume->IsPoly())
 918   {
 919     myPolyIndices.resize( myFaceNbNodes + 1 );
 920     myFaceNodeIndices = & myPolyIndices[0];
 921     for ( int i = 0; i <= myFaceNbNodes; ++i )
 922       myFaceNodeIndices[i] = myVolume->GetNodeIndex( myFaceNodes[i] );
 923   }
 924   return myFaceNodeIndices;
 925 }
 926
 927 //=======================================================================
 928 //function : GetFaceNodes
 929 //purpose  : Return a set of face nodes.
 930 //=======================================================================
 931
 932 bool SMDS_VolumeTool::GetFaceNodes (int                        faceIndex,
 933                                     set<const SMDS_MeshNode*>& theFaceNodes )
 934 {
 935   if ( !setFace( faceIndex ))
 936     return false;
 937
 938   theFaceNodes.clear();
 939   int iNode, nbNode = myFaceNbNodes;
 940   for ( iNode = 0; iNode < nbNode; iNode++ )
 941     theFaceNodes.insert( myFaceNodes[ iNode ]);
 942
 943   return true;
 944 }
 945
 946 //=======================================================================
 947 //function : IsFaceExternal
 948 //purpose  : Check normal orientation of a returned face
 949 //=======================================================================
 950
 951 bool SMDS_VolumeTool::IsFaceExternal( int faceIndex )
 952 {
 953   if ( myExternalFaces || !myVolume )
 954     return true;
 955
 956   if (myVolume->IsPoly()) {
 957     XYZ aNormal, baryCenter, p0 (myPolyedre->GetFaceNode(faceIndex + 1, 1));
 958     GetFaceNormal(faceIndex, aNormal.x, aNormal.y, aNormal.z);
 959     GetBaryCenter(baryCenter.x, baryCenter.y, baryCenter.z);
 960     XYZ insideVec (baryCenter - p0);
 961     if (insideVec.Dot(aNormal) > 0)
 962       return false;
 963     return true;
 964   }
 965
 966   switch ( myVolumeNbNodes ) {
 967   case 4:
 968   case 5:
 969   case 10:
 970   case 13:
 971     // only the bottom of a reversed tetrahedron can be internal
 972     return ( myVolForward || faceIndex != 0 );
 973   case 6:
 974   case 15:
 975     // in a forward pentahedron, the top is internal, in a reversed one - bottom
 976     return ( myVolForward ? faceIndex != 1 : faceIndex != 0 );
 977   case 8:
 978   case 20: {
 979     // in a forward hexahedron, even face normal is external, odd - internal
 980     bool odd = faceIndex % 2;
 981     return ( myVolForward ? !odd : odd );
 982   }
 983   default:;
 984   }
 985   return false;
 986 }
 987
 988 //=======================================================================
 989 //function : GetFaceNormal
 990 //purpose  : Return a normal to a face
 991 //=======================================================================
 992
 993 bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z)
 994 {
 995   if ( !setFace( faceIndex ))
 996     return false;
 997
 998   XYZ p1 ( myFaceNodes[0] );
 999   XYZ p2 ( myFaceNodes[1] );
1000   XYZ p3 ( myFaceNodes[2] );
1001   XYZ aVec12( p2 - p1 );
1002   XYZ aVec13( p3 - p1 );
1003   XYZ cross = aVec12.Crossed( aVec13 );
1004
1005   //if ( myFaceNbNodes == 4 ) {
1006   if ( myFaceNbNodes >3 ) {
1007     XYZ p4 ( myFaceNodes[3] );
1008     XYZ aVec14( p4 - p1 );
1009     XYZ cross2 = aVec13.Crossed( aVec14 );
1010     cross.x += cross2.x;
1011     cross.y += cross2.y;
1012     cross.z += cross2.z;
1013   }
1014
1015   double size = cross.Magnitude();
1016   if ( size <= DBL_MIN )
1017     return false;
1018
1019   X = cross.x / size;
1020   Y = cross.y / size;
1021   Z = cross.z / size;
1022
1023   return true;
1024 }
1025
1026 //=======================================================================
1027 //function : GetFaceArea
1028 //purpose  : Return face area
1029 //=======================================================================
1030
1031 double SMDS_VolumeTool::GetFaceArea( int faceIndex )
1032 {
1033   if (myVolume->IsPoly()) {
1034     MESSAGE("Warning: attempt to obtain area of a face of polyhedral volume");
1035     return 0;
1036   }
1037
1038   if ( !setFace( faceIndex ))
1039     return 0;
1040
1041   XYZ p1 ( myFaceNodes[0] );
1042   XYZ p2 ( myFaceNodes[1] );
1043   XYZ p3 ( myFaceNodes[2] );
1044   XYZ aVec12( p2 - p1 );
1045   XYZ aVec13( p3 - p1 );
1046   double area = aVec12.Crossed( aVec13 ).Magnitude() * 0.5;
1047
1048   if ( myFaceNbNodes == 4 ) {
1049     XYZ p4 ( myFaceNodes[3] );
1050     XYZ aVec14( p4 - p1 );
1051     area += aVec14.Crossed( aVec13 ).Magnitude() * 0.5;
1052   }
1053   return area;
1054 }
1055
1056 //=======================================================================
1057 //function : GetOppFaceIndex
1058 //purpose  : Return index of the opposite face if it exists, else -1.
1059 //=======================================================================
1060
1061 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
1062 {
1063   int ind = -1;
1064   if (myVolume->IsPoly()) {
1065     MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume");
1066     return ind;
1067   }
1068
1069   if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
1070     switch ( myVolumeNbNodes ) {
1071     case 6:
1072       if ( faceIndex == 0 || faceIndex == 1 )
1073         ind = 1 - faceIndex;
1074         break;
1075     case 8:
1076       ind = faceIndex + ( faceIndex % 2 ? -1 : 1 );
1077       break;
1078     default:;
1079     }
1080   }
1081   return ind;
1082 }
1083
1084 //=======================================================================
1085 //function : IsLinked
1086 //purpose  : return true if theNode1 is linked with theNode2
1087 //=======================================================================
1088
1089 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
1090                                 const SMDS_MeshNode* theNode2) const
1091 {
1092   if ( !myVolume )
1093     return false;
1094
1095   if (myVolume->IsPoly()) {
1096     if (!myPolyedre) {
1097       MESSAGE("Warning: bad volumic element");
1098       return false;
1099     }
1100     bool isLinked = false;
1101     int iface;
1102     for (iface = 1; iface <= myNbFaces && !isLinked; iface++) {
1103       int inode, nbFaceNodes = myPolyedre->NbFaceNodes(iface);
1104
1105       for (inode = 1; inode <= nbFaceNodes && !isLinked; inode++) {
1106         const SMDS_MeshNode* curNode = myPolyedre->GetFaceNode(iface, inode);
1107
1108         if (curNode == theNode1 || curNode == theNode2) {
1109           int inextnode = (inode == nbFaceNodes) ? 1 : inode + 1;
1110           const SMDS_MeshNode* nextNode = myPolyedre->GetFaceNode(iface, inextnode);
1111
1112           if ((curNode == theNode1 && nextNode == theNode2) ||
1113               (curNode == theNode2 && nextNode == theNode1)) {
1114             isLinked = true;
1115           }
1116         }
1117       }
1118     }
1119     return isLinked;
1120   }
1121
1122   // find nodes indices
1123   int i1 = -1, i2 = -1;
1124   for ( int i = 0; i < myVolumeNbNodes; i++ ) {
1125     if ( myVolumeNodes[ i ] == theNode1 )
1126       i1 = i;
1127     else if ( myVolumeNodes[ i ] == theNode2 )
1128       i2 = i;
1129   }
1130   return IsLinked( i1, i2 );
1131 }
1132
1133 //=======================================================================
1134 //function : IsLinked
1135 //purpose  : return true if the node with theNode1Index is linked
1136 //           with the node with theNode2Index
1137 //=======================================================================
1138
1139 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
1140                                 const int theNode2Index) const
1141 {
1142   if ( myVolume->IsPoly() ) {
1143     return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
1144   }
1145
1146   int minInd = min( theNode1Index, theNode2Index );
1147   int maxInd = max( theNode1Index, theNode2Index );
1148
1149   if ( minInd < 0 || maxInd > myVolumeNbNodes - 1 || maxInd == minInd )
1150     return false;
1151
1152   switch ( myVolumeNbNodes ) {
1153   case 4:
1154     return true;
1155   case 5:
1156     if ( maxInd == 4 )
1157       return true;
1158     switch ( maxInd - minInd ) {
1159     case 1:
1160     case 3: return true;
1161     default:;
1162     }
1163     break;
1164   case 6:
1165     switch ( maxInd - minInd ) {
1166     case 1: return minInd != 2;
1167     case 2: return minInd == 0 || minInd == 3;
1168     case 3: return true;
1169     default:;
1170     }
1171     break;
1172   case 8:
1173     switch ( maxInd - minInd ) {
1174     case 1: return minInd != 3;
1175     case 3: return minInd == 0 || minInd == 4;
1176     case 4: return true;
1177     default:;
1178     }
1179     break;
1180   case 10:
1181     {
1182       switch ( minInd ) {
1183       case 0: if( maxInd==4 ||  maxInd==6 ||  maxInd==7 ) return true;
1184       case 1: if( maxInd==4 ||  maxInd==5 ||  maxInd==8 ) return true;
1185       case 2: if( maxInd==5 ||  maxInd==6 ||  maxInd==9 ) return true;
1186       case 3: if( maxInd==7 ||  maxInd==8 ||  maxInd==9 ) return true;
1187       default:;
1188       }
1189       break;
1190     }
1191   case 13:
1192     {
1193       switch ( minInd ) {
1194       case 0: if( maxInd==5 ||  maxInd==8 ||  maxInd==9 ) return true;
1195       case 1: if( maxInd==5 ||  maxInd==6 ||  maxInd==10 ) return true;
1196       case 2: if( maxInd==6 ||  maxInd==7 ||  maxInd==11 ) return true;
1197       case 3: if( maxInd==7 ||  maxInd==8 ||  maxInd==12 ) return true;
1198       case 4: if( maxInd==9 ||  maxInd==10 ||  maxInd==11 ||  maxInd==12 ) return true;
1199       default:;
1200       }
1201       break;
1202     }
1203   case 15:
1204     {
1205       switch ( minInd ) {
1206       case 0: if( maxInd==6 ||  maxInd==8 ||  maxInd==12 ) return true;
1207       case 1: if( maxInd==6 ||  maxInd==7 ||  maxInd==13 ) return true;
1208       case 2: if( maxInd==7 ||  maxInd==8 ||  maxInd==14 ) return true;
1209       case 3: if( maxInd==9 ||  maxInd==11 ||  maxInd==12 ) return true;
1210       case 4: if( maxInd==9 ||  maxInd==10 ||  maxInd==13 ) return true;
1211       case 5: if( maxInd==10 ||  maxInd==11 ||  maxInd==14 ) return true;
1212       default:;
1213       }
1214       break;
1215     }
1216   case 20:
1217     {
1218       switch ( minInd ) {
1219       case 0: if( maxInd==8 ||  maxInd==11 ||  maxInd==16 ) return true;
1220       case 1: if( maxInd==8 ||  maxInd==9 ||  maxInd==17 ) return true;
1221       case 2: if( maxInd==9 ||  maxInd==10 ||  maxInd==18 ) return true;
1222       case 3: if( maxInd==10 ||  maxInd==11 ||  maxInd==19 ) return true;
1223       case 4: if( maxInd==12 ||  maxInd==15 ||  maxInd==16 ) return true;
1224       case 5: if( maxInd==12 ||  maxInd==13 ||  maxInd==17 ) return true;
1225       case 6: if( maxInd==13 ||  maxInd==14 ||  maxInd==18 ) return true;
1226       case 7: if( maxInd==14 ||  maxInd==15 ||  maxInd==19 ) return true;
1227       default:;
1228       }
1229       break;
1230     }
1231   default:;
1232   }
1233   return false;
1234 }
1235
1236 //=======================================================================
1237 //function : GetNodeIndex
1238 //purpose  : Return an index of theNode
1239 //=======================================================================
1240
1241 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
1242 {
1243   if ( myVolume ) {
1244     for ( int i = 0; i < myVolumeNbNodes; i++ ) {
1245       if ( myVolumeNodes[ i ] == theNode )
1246         return i;
1247     }
1248   }
1249   return -1;
1250 }
1251
1252 //================================================================================
1253 /*!
1254  * \brief Fill vector with boundary faces existing in the mesh
1255   * \param faces - vector of found nodes
1256   * \retval int - nb of found faces
1257  */
1258 //================================================================================
1259
1260 int SMDS_VolumeTool::GetAllExistingFaces(vector<const SMDS_MeshElement*> & faces)
1261 {
1262   faces.clear();
1263   faces.reserve( NbFaces() );
1264   for ( int iF = 0; iF < NbFaces(); ++iF ) {
1265     const SMDS_MeshFace* face = 0;
1266     const SMDS_MeshNode** nodes = GetFaceNodes( iF );
1267     switch ( NbFaceNodes( iF )) {
1268     case 3:
1269       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2] ); break;
1270     case 4:
1271       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
1272     case 6:
1273       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2],
1274                                   nodes[3], nodes[4], nodes[5]); break;
1275     case 8:
1276       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3],
1277                                   nodes[4], nodes[5], nodes[6], nodes[7]); break;
1278     }
1279     if ( face )
1280       faces.push_back( face );
1281   }
1282   return faces.size();
1283 }
1284
1285
1286 //================================================================================
1287 /*!
1288  * \brief Fill vector with boundary edges existing in the mesh
1289   * \param edges - vector of found edges
1290   * \retval int - nb of found faces
1291  */
1292 //================================================================================
1293
1294 int SMDS_VolumeTool::GetAllExistingEdges(vector<const SMDS_MeshElement*> & edges) const
1295 {
1296   edges.clear();
1297   edges.reserve( myVolumeNbNodes * 2 );
1298   for ( int i = 0; i < myVolumeNbNodes; ++i ) {
1299     for ( int j = i + 1; j < myVolumeNbNodes; ++j ) {
1300       if ( IsLinked( i, j )) {
1301         const SMDS_MeshElement* edge =
1302           SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] );
1303         if ( edge )
1304           edges.push_back( edge );
1305       }
1306     }
1307   }
1308   return edges.size();
1309 }
1310
1311 //=======================================================================
1312 //function : IsFreeFace
1313 //purpose  : check that only one volume is build on the face nodes
1314 //=======================================================================
1315
1316 bool SMDS_VolumeTool::IsFreeFace( int faceIndex )
1317 {
1318   const int free = true;
1319
1320   if (!setFace( faceIndex ))
1321     return !free;
1322
1323   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1324   int nbFaceNodes = myFaceNbNodes;
1325
1326   // evaluate nb of face nodes shared by other volume
1327   int maxNbShared = -1;
1328   typedef map< const SMDS_MeshElement*, int > TElemIntMap;
1329   TElemIntMap volNbShared;
1330   TElemIntMap::iterator vNbIt;
1331   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1332     const SMDS_MeshNode* n = nodes[ iNode ];
1333     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator();
1334     while ( eIt->more() ) {
1335       const SMDS_MeshElement* elem = eIt->next();
1336       if ( elem != myVolume && elem->GetType() == SMDSAbs_Volume ) {
1337         int nbShared = 1;
1338         vNbIt = volNbShared.find( elem );
1339         if ( vNbIt == volNbShared.end() ) {
1340           volNbShared.insert ( TElemIntMap::value_type( elem, nbShared ));
1341         }
1342         else {
1343           nbShared = ++(*vNbIt).second;
1344         }
1345         if ( nbShared > maxNbShared )
1346           maxNbShared = nbShared;
1347       }
1348     }
1349   }
1350   if ( maxNbShared < 3 )
1351     return free; // is free
1352
1353   // find volumes laying on the opposite side of the face
1354   // and sharing all nodes
1355   XYZ intNormal; // internal normal
1356   GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
1357   if ( IsFaceExternal( faceIndex ))
1358     intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
1359   XYZ p0 ( nodes[0] ), baryCenter;
1360   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
1361     int nbShared = (*vNbIt).second;
1362     if ( nbShared >= 3 ) {
1363       SMDS_VolumeTool volume( (*vNbIt).first );
1364       volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
1365       XYZ intNormal2( baryCenter - p0 );
1366       if ( intNormal.Dot( intNormal2 ) < 0 )
1367         continue; // opposite side
1368     }
1369     // remove a volume from volNbShared map
1370     volNbShared.erase( vNbIt-- );
1371   }
1372
1373   // here volNbShared contains only volumes laying on the
1374   // opposite side of the face
1375   if ( volNbShared.empty() ) {
1376     return free; // is free
1377   }
1378
1379   // check if the whole area of a face is shared
1380   bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
1381   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
1382     SMDS_VolumeTool volume( (*vNbIt).first );
1383     bool prevLinkShared = false;
1384     int nbSharedLinks = 0;
1385     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1386       bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
1387       if ( linkShared )
1388         nbSharedLinks++;
1389       if ( linkShared && prevLinkShared &&
1390           volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
1391         isShared[ iNode ] = true;
1392       prevLinkShared = linkShared;
1393     }
1394     if ( nbSharedLinks == nbFaceNodes )
1395       return !free; // is not free
1396     if ( nbFaceNodes == 4 ) {
1397       // check traingle parts 1 & 3
1398       if ( isShared[1] && isShared[3] )
1399         return !free; // is not free
1400       // check triangle parts 0 & 2;
1401       // 0 part could not be checked in the loop; check it here
1402       if ( isShared[2] && prevLinkShared &&
1403           volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
1404           volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
1405         return !free; // is not free
1406     }
1407   }
1408   return free;
1409 }
1410
1411 //=======================================================================
1412 //function : GetFaceIndex
1413 //purpose  : Return index of a face formed by theFaceNodes
1414 //=======================================================================
1415
1416 int SMDS_VolumeTool::GetFaceIndex( const set<const SMDS_MeshNode*>& theFaceNodes )
1417 {
1418   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
1419     const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
1420     int nbFaceNodes = NbFaceNodes( iFace );
1421     set<const SMDS_MeshNode*> nodeSet;
1422     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1423       nodeSet.insert( nodes[ iNode ] );
1424     if ( theFaceNodes == nodeSet )
1425       return iFace;
1426   }
1427   return -1;
1428 }
1429
1430 //=======================================================================
1431 //function : GetFaceIndex
1432 //purpose  : Return index of a face formed by theFaceNodes
1433 //=======================================================================
1434
1435 /*int SMDS_VolumeTool::GetFaceIndex( const set<int>& theFaceNodesIndices )
1436 {
1437   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
1438     const int* nodes = GetFaceNodesIndices( iFace );
1439     int nbFaceNodes = NbFaceNodes( iFace );
1440     set<int> nodeSet;
1441     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1442       nodeSet.insert( nodes[ iNode ] );
1443     if ( theFaceNodesIndices == nodeSet )
1444       return iFace;
1445   }
1446   return -1;
1447 }*/
1448
1449 //=======================================================================
1450 //function : setFace
1451 //purpose  :
1452 //=======================================================================
1453
1454 bool SMDS_VolumeTool::setFace( int faceIndex )
1455 {
1456   if ( !myVolume )
1457     return false;
1458
1459   if ( myCurFace == faceIndex )
1460     return true;
1461
1462   myCurFace = -1;
1463
1464   if ( faceIndex < 0 || faceIndex >= NbFaces() )
1465     return false;
1466
1467   if (myFaceNodes != NULL) {
1468     delete [] myFaceNodes;
1469     myFaceNodes = NULL;
1470   }
1471
1472   if (myVolume->IsPoly()) {
1473     if (!myPolyedre) {
1474       MESSAGE("Warning: bad volumic element");
1475       return false;
1476     }
1477
1478     // check orientation
1479     bool isGoodOri = true;
1480     if (myExternalFaces)
1481       isGoodOri = IsFaceExternal( faceIndex );
1482
1483     // set face nodes
1484     int iNode;
1485     myFaceNbNodes = myPolyedre->NbFaceNodes(faceIndex + 1);
1486     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
1487     if (isGoodOri) {
1488       for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
1489         myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, iNode + 1);
1490     } else {
1491       for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
1492         myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, myFaceNbNodes - iNode);
1493     }
1494     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ]; // last = first
1495
1496   }
1497   else {
1498     // choose face node indices
1499     switch ( myVolumeNbNodes ) {
1500     case 4:
1501       myFaceNbNodes = Tetra_nbN[ faceIndex ];
1502       if ( myExternalFaces )
1503         myFaceNodeIndices = myVolForward ? Tetra_F[ faceIndex ] : Tetra_RE[ faceIndex ];
1504       else
1505         myFaceNodeIndices = myVolForward ? Tetra_F[ faceIndex ] : Tetra_R[ faceIndex ];
1506       break;
1507     case 5:
1508       myFaceNbNodes = Pyramid_nbN[ faceIndex ];
1509       if ( myExternalFaces )
1510         myFaceNodeIndices = myVolForward ? Pyramid_F[ faceIndex ] : Pyramid_RE[ faceIndex ];
1511       else
1512         myFaceNodeIndices = myVolForward ? Pyramid_F[ faceIndex ] : Pyramid_R[ faceIndex ];
1513       break;
1514     case 6:
1515       myFaceNbNodes = Penta_nbN[ faceIndex ];
1516       if ( myExternalFaces )
1517         myFaceNodeIndices = myVolForward ? Penta_FE[ faceIndex ] : Penta_RE[ faceIndex ];
1518       else
1519         myFaceNodeIndices = myVolForward ? Penta_F[ faceIndex ] : Penta_R[ faceIndex ];
1520       break;
1521     case 8:
1522       myFaceNbNodes = Hexa_nbN[ faceIndex ];
1523       if ( myExternalFaces )
1524         myFaceNodeIndices = myVolForward ? Hexa_FE[ faceIndex ] : Hexa_RE[ faceIndex ];
1525       else
1526         myFaceNodeIndices = Hexa_F[ faceIndex ];
1527       break;
1528     case 10:
1529       myFaceNbNodes = QuadTetra_nbN[ faceIndex ];
1530       if ( myExternalFaces )
1531         myFaceNodeIndices = myVolForward ? QuadTetra_F[ faceIndex ] : QuadTetra_RE[ faceIndex ];
1532       else
1533         myFaceNodeIndices = myVolForward ? QuadTetra_F[ faceIndex ] : QuadTetra_R[ faceIndex ];
1534       break;
1535     case 13:
1536       myFaceNbNodes = QuadPyram_nbN[ faceIndex ];
1537       if ( myExternalFaces )
1538         myFaceNodeIndices = myVolForward ? QuadPyram_F[ faceIndex ] : QuadPyram_RE[ faceIndex ];
1539       else
1540         myFaceNodeIndices = myVolForward ? QuadPyram_F[ faceIndex ] : QuadPyram_R[ faceIndex ];
1541       break;
1542     case 15:
1543       myFaceNbNodes = QuadPenta_nbN[ faceIndex ];
1544       if ( myExternalFaces )
1545         myFaceNodeIndices = myVolForward ? QuadPenta_FE[ faceIndex ] : QuadPenta_RE[ faceIndex ];
1546       else
1547         myFaceNodeIndices = myVolForward ? QuadPenta_F[ faceIndex ] : QuadPenta_R[ faceIndex ];
1548       break;
1549     case 20:
1550       myFaceNbNodes = QuadHexa_nbN[ faceIndex ];
1551       if ( myExternalFaces )
1552         myFaceNodeIndices = myVolForward ? QuadHexa_FE[ faceIndex ] : QuadHexa_RE[ faceIndex ];
1553       else
1554         myFaceNodeIndices = QuadHexa_F[ faceIndex ];
1555       break;
1556     default:
1557       return false;
1558     }
1559
1560     // set face nodes
1561     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
1562     for ( int iNode = 0; iNode < myFaceNbNodes; iNode++ )
1563       myFaceNodes[ iNode ] = myVolumeNodes[ myFaceNodeIndices[ iNode ]];
1564     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ];
1565   }
1566
1567   myCurFace = faceIndex;
1568
1569   return true;
1570 }
1571
1572 //=======================================================================
1573 //function : GetType
1574 //purpose  : return VolumeType by nb of nodes in a volume
1575 //=======================================================================
1576
1577 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
1578 {
1579   switch ( nbNodes ) {
1580   case 4: return TETRA;
1581   case 5: return PYRAM;
1582   case 6: return PENTA;
1583   case 8: return HEXA;
1584   case 10: return QUAD_TETRA;
1585   case 13: return QUAD_PYRAM;
1586   case 15: return QUAD_PENTA;
1587   case 20: return QUAD_HEXA;
1588   default:return UNKNOWN;
1589   }
1590 }
1591
1592 //=======================================================================
1593 //function : NbFaces
1594 //purpose  : return nb of faces by volume type
1595 //=======================================================================
1596
1597 int SMDS_VolumeTool::NbFaces( VolumeType type )
1598 {
1599   switch ( type ) {
1600   case TETRA     :
1601   case QUAD_TETRA: return 4;
1602   case PYRAM     :
1603   case QUAD_PYRAM: return 5;
1604   case PENTA     :
1605   case QUAD_PENTA: return 5;
1606   case HEXA      :
1607   case QUAD_HEXA : return 6;
1608   default:    return 0;
1609   }
1610 }
1611
1612 //================================================================================
1613 /*!
1614  * \brief Useful to know nb of corner nodes of a quadratic volume
1615   * \param type - volume type
1616   * \retval int - nb of corner nodes
1617  */
1618 //================================================================================
1619
1620 int SMDS_VolumeTool::NbCornerNodes(VolumeType type)
1621 {
1622   switch ( type ) {
1623   case TETRA     :
1624   case QUAD_TETRA: return 4;
1625   case PYRAM     :
1626   case QUAD_PYRAM: return 5;
1627   case PENTA     :
1628   case QUAD_PENTA: return 6;
1629   case HEXA      :
1630   case QUAD_HEXA : return 8;
1631   default:    return 0;
1632   }
1633   return 0;
1634 }
1635   //
1636
1637 //=======================================================================
1638 //function : GetFaceNodesIndices
1639 //purpose  : Return the array of face nodes indices
1640 //           To comfort link iteration, the array
1641 //           length == NbFaceNodes( faceIndex ) + 1 and
1642 //           the last node index == the first one.
1643 //=======================================================================
1644
1645 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
1646                                                 int        faceIndex,
1647                                                 bool       external)
1648 {
1649   switch ( type ) {
1650   case TETRA: return Tetra_F[ faceIndex ];
1651   case PYRAM: return Pyramid_F[ faceIndex ];
1652   case PENTA: return external ? Penta_FE[ faceIndex ] : Penta_F[ faceIndex ];
1653   case HEXA:  return external ? Hexa_FE[ faceIndex ] : Hexa_F[ faceIndex ];
1654   case QUAD_TETRA: return QuadTetra_F[ faceIndex ];
1655   case QUAD_PYRAM: return QuadPyram_F[ faceIndex ];
1656   case QUAD_PENTA: return external ? QuadPenta_FE[ faceIndex ] : QuadPenta_F[ faceIndex ];
1657   case QUAD_HEXA:  return external ? QuadHexa_FE[ faceIndex ] : QuadHexa_F[ faceIndex ];
1658   default:;
1659   }
1660   return 0;
1661 }
1662
1663 //=======================================================================
1664 //function : NbFaceNodes
1665 //purpose  : Return number of nodes in the array of face nodes
1666 //=======================================================================
1667
1668 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
1669                                  int        faceIndex )
1670 {
1671   switch ( type ) {
1672   case TETRA: return Tetra_nbN[ faceIndex ];
1673   case PYRAM: return Pyramid_nbN[ faceIndex ];
1674   case PENTA: return Penta_nbN[ faceIndex ];
1675   case HEXA:  return Hexa_nbN[ faceIndex ];
1676   case QUAD_TETRA: return QuadTetra_nbN[ faceIndex ];
1677   case QUAD_PYRAM: return QuadPyram_nbN[ faceIndex ];
1678   case QUAD_PENTA: return QuadPenta_nbN[ faceIndex ];
1679   case QUAD_HEXA:  return QuadHexa_nbN[ faceIndex ];
1680   default:;
1681   }
1682   return 0;
1683 }
1684
1685 //=======================================================================
1686 //function : Get
1687 //purpose  : return element
1688 //=======================================================================
1689
1690 const SMDS_MeshVolume* SMDS_VolumeTool::Element() const
1691 {
1692   return static_cast<const SMDS_MeshVolume*>( myVolume );
1693 }
1694
1695 //=======================================================================
1696 //function : ID
1697 //purpose  : return element ID
1698 //=======================================================================
1699
1700 int SMDS_VolumeTool::ID() const
1701 {
1702   return myVolume ? myVolume->GetID() : 0;
1703 }