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