src/SMDS/SMDS_VolumeTool.cxx

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