src/SMDS/SMDS_VolumeTool.cxx

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