Salome HOME
This commit was generated by cvs2git to create tag 'T2_2_5a'.
[modules/smesh.git] / src / SMDS / SMDS_VolumeTool.cxx
1 // File      : SMDS_VolumeTool.cxx
2 // Created   : Tue Jul 13 12:22:13 2004
3 // Author    : Edward AGAPOV (eap)
4 // Copyright : Open CASCADE
5
6 #ifdef _MSC_VER
7 #pragma warning(disable:4786)
8 #endif
9
10 #include "SMDS_VolumeTool.hxx"
11
12 #include "SMDS_MeshElement.hxx"
13 #include "SMDS_MeshNode.hxx"
14 #include <map>
15 #include <float.h>
16 #include <math.h>
17
18 using namespace std;
19
20 // ======================================================
21 // Node indices in faces depending on volume orientation
22 // making most faces normals external
23 // ======================================================
24
25 /*
26 //           N3
27 //           +
28 //          /|\
29 //         / | \
30 //        /  |  \
31 //    N0 +---|---+ N1                TETRAHEDRON
32 //       \   |   /
33 //        \  |  /
34 //         \ | /
35 //          \|/
36 //           +
37 //           N2
38 */
39 static int Tetra_F [4][4] = { // FORWARD == EXTERNAL
40   { 0, 1, 2, 0 },              // All faces have external normals
41   { 0, 3, 1, 0 },
42   { 1, 3, 2, 1 },
43   { 0, 2, 3, 0 }}; 
44 static int Tetra_R [4][4] = { // REVERSED
45   { 0, 1, 2, 0 },             // All faces but a bottom have external normals
46   { 0, 1, 3, 0 },
47   { 1, 2, 3, 1 },
48   { 0, 3, 2, 0 }};
49 static int Tetra_RE [4][4] = { // REVERSED -> FORWARD (EXTERNAL)
50   { 0, 2, 1, 0 },              // All faces have external normals
51   { 0, 1, 3, 0 },
52   { 1, 2, 3, 1 },
53   { 0, 3, 2, 0 }};
54 static int Tetra_nbN [] = { 3, 3, 3, 3 };
55
56 //
57 //     PYRAMID
58 //
59 static int Pyramid_F [5][5] = { // FORWARD == EXTERNAL
60   { 0, 1, 2, 3, 0 },            // All faces have external normals
61   { 0, 4, 1, 0, 4 },
62   { 1, 4, 2, 1, 4 },
63   { 2, 4, 3, 2, 4 },
64   { 3, 4, 0, 3, 4 }}; 
65 static int Pyramid_R [5][5] = { // REVERSED
66   { 0, 1, 2, 3, 0 },            // All faces but a bottom have external normals
67   { 0, 1, 4, 0, 4 },
68   { 1, 2, 4, 1, 4 },
69   { 2, 3, 4, 2, 4 },
70   { 3, 0, 4, 3, 4 }}; 
71 static int Pyramid_RE [5][5] = { // REVERSED -> FORWARD (EXTERNAL)
72   { 0, 3, 2, 1, 0 },             // All faces but a bottom have external normals
73   { 0, 1, 4, 0, 4 },
74   { 1, 2, 4, 1, 4 },
75   { 2, 3, 4, 2, 4 },
76   { 3, 0, 4, 3, 4 }}; 
77 static int Pyramid_nbN [] = { 4, 3, 3, 3, 3 };
78
79 /*   
80 //            + N4
81 //           /|\
82 //          / | \
83 //         /  |  \
84 //        /   |   \
85 //    N3 +---------+ N5
86 //       |    |    |
87 //       |    + N1 |
88 //       |   / \   |                PENTAHEDRON
89 //       |  /   \  |
90 //       | /     \ |
91 //       |/       \|
92 //    N0 +---------+ N2
93 */
94 static int Penta_F [5][5] = { // FORWARD
95   { 0, 1, 2, 0, 0 },          // Top face has an internal normal, other - external
96   { 3, 4, 5, 3, 3 },          // 0 is bottom, 1 is top face
97   { 0, 2, 5, 3, 0 },
98   { 1, 4, 5, 2, 1 },
99   { 0, 3, 4, 1, 0 }}; 
100 static int Penta_R [5][5] = { // REVERSED
101   { 0, 1, 2, 0, 0 },          // Bottom face has an internal normal, other - external
102   { 3, 4, 5, 3, 3 },          // 0 is bottom, 1 is top face
103   { 0, 3, 5, 2, 0 },
104   { 1, 2, 5, 4, 1 },
105   { 0, 1, 4, 3, 0 }}; 
106 static int Penta_FE [5][5] = { // FORWARD -> EXTERNAL
107   { 0, 1, 2, 0, 0 },
108   { 3, 5, 4, 3, 3 },
109   { 0, 2, 5, 3, 0 },
110   { 1, 4, 5, 2, 1 },
111   { 0, 3, 4, 1, 0 }}; 
112 static int Penta_RE [5][5] = { // REVERSED -> EXTERNAL
113   { 0, 2, 1, 0, 0 },
114   { 3, 4, 5, 3, 3 },
115   { 0, 3, 5, 2, 0 },
116   { 1, 2, 5, 4, 1 },
117   { 0, 1, 4, 3, 0 }}; 
118 static int Penta_nbN [] = { 3, 3, 4, 4, 4 };
119
120 /*
121 //         N5+----------+N6
122 //          /|         /|
123 //         / |        / |
124 //        /  |       /  |
125 //     N4+----------+N7 |
126 //       |   |      |   |           HEXAHEDRON
127 //       |   |      |   |
128 //       |   |      |   |
129 //       | N1+------|---+N2
130 //       |  /       |  /
131 //       | /        | /
132 //       |/         |/
133 //     N0+----------+N3
134 */
135 static int Hexa_F [6][5] = { // FORWARD
136   { 0, 1, 2, 3, 0 },         // opposite faces are neighbouring,
137   { 4, 5, 6, 7, 4 },         // odd face(1,3,5) normal is internal, even(0,2,4) - external
138   { 1, 0, 4, 5, 1 },         // same index nodes of opposite faces are linked
139   { 2, 3, 7, 6, 2 }, 
140   { 0, 3, 7, 4, 0 }, 
141   { 1, 2, 6, 5, 1 }};
142 // static int Hexa_R [6][5] = { // REVERSED
143 //   { 0, 3, 2, 1, 0 },         // opposite faces are neighbouring,
144 //   { 4, 7, 6, 5, 4 },         // odd face(1,3,5) normal is external, even(0,2,4) - internal
145 //   { 1, 5, 4, 0, 1 },         // same index nodes of opposite faces are linked
146 //   { 2, 6, 7, 3, 2 }, 
147 //   { 0, 4, 7, 3, 0 }, 
148 //   { 1, 5, 6, 2, 1 }};
149 static int Hexa_FE [6][5] = { // FORWARD -> EXTERNAL
150   { 0, 1, 2, 3, 0 } ,         // opposite faces are neighbouring,
151   { 4, 7, 6, 5, 4 },          // all face normals are external,
152   { 0, 4, 5, 1, 0 },          // links in opposite faces: 0-0, 1-3, 2-2, 3-1
153   { 3, 2, 6, 7, 3 }, 
154   { 0, 3, 7, 4, 0 },
155   { 1, 5, 6, 2, 1 }};
156 static int Hexa_RE [6][5] = { // REVERSED -> EXTERNAL
157   { 0, 3, 2, 1, 0 },          // opposite faces are neighbouring,
158   { 4, 5, 6, 7, 4 },          // all face normals are external,
159   { 0, 1, 5, 4, 0 },          // links in opposite faces: 0-0, 1-3, 2-2, 3-1
160   { 3, 7, 6, 2, 3 }, 
161   { 0, 4, 7, 3, 0 },
162   { 1, 2, 6, 5, 1 }};
163 static int Hexa_nbN [] = { 4, 4, 4, 4, 4, 4 };
164
165 // ========================================================
166 // to perform some calculations without linkage to CASCADE
167 // ========================================================
168 struct XYZ {
169   double x;
170   double y;
171   double z;
172   XYZ()                               { x = 0; y = 0; z = 0; }
173   XYZ( double X, double Y, double Z ) { x = X; y = Y; z = Z; }
174   XYZ( const XYZ& other )             { x = other.x; y = other.y; z = other.z; }
175   XYZ( const SMDS_MeshNode* n )       { x = n->X(); y = n->Y(); z = n->Z(); }
176   XYZ operator-( const XYZ& other );
177   XYZ Crossed( const XYZ& other );
178   double Dot( const XYZ& other );
179   double Magnitude();
180 };
181 XYZ XYZ::operator-( const XYZ& Right ) {
182   return XYZ(x - Right.x, y - Right.y, z - Right.z);
183 }
184 XYZ XYZ::Crossed( const XYZ& Right ) {
185   return XYZ (y * Right.z - z * Right.y,
186               z * Right.x - x * Right.z,
187               x * Right.y - y * Right.x);
188 }
189 double XYZ::Dot( const XYZ& Other ) {
190   return(x * Other.x + y * Other.y + z * Other.z);
191 }
192 double XYZ::Magnitude() {
193   return sqrt (x * x + y * y + z * z);
194 }
195
196 //=======================================================================
197 //function : SMDS_VolumeTool
198 //purpose  : 
199 //=======================================================================
200
201 SMDS_VolumeTool::SMDS_VolumeTool ()
202      : myVolume( 0 ),
203        myVolForward( true ),
204        myNbFaces( 0 ),
205        myVolumeNbNodes( 0 ),
206        myExternalFaces( false )
207 {
208 }
209 //=======================================================================
210 //function : SMDS_VolumeTool
211 //purpose  : 
212 //=======================================================================
213
214 SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume)
215      : myExternalFaces( false )
216 {
217   Set( theVolume );
218 }
219
220 //=======================================================================
221 //function : SMDS_VolumeTool
222 //purpose  : 
223 //=======================================================================
224
225 SMDS_VolumeTool::~SMDS_VolumeTool()
226 {
227 }
228
229 //=======================================================================
230 //function : SetVolume
231 //purpose  : Set volume to iterate on
232 //=======================================================================
233
234 bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume)
235 {
236   myVolume = 0;
237   myVolForward = true;
238   myCurFace = -1;
239   myVolumeNbNodes = 0;
240   myNbFaces = 0;
241   if ( theVolume && theVolume->GetType() == SMDSAbs_Volume )
242   {
243     myVolumeNbNodes = theVolume->NbNodes();
244     switch ( myVolumeNbNodes ) {
245     case 4:
246     case 5:
247     case 6:
248     case 8:
249       {
250       myVolume = theVolume;
251       myNbFaces = theVolume->NbFaces();
252
253       // set volume nodes
254       int iNode = 0;
255       SMDS_ElemIteratorPtr nodeIt = myVolume->nodesIterator();
256       while ( nodeIt->more() )
257         myVolumeNodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
258
259       // nb nodes in each face
260       if ( myVolumeNbNodes == 4 )
261         myFaceNbNodes = Tetra_nbN;
262       else if ( myVolumeNbNodes == 5 )
263         myFaceNbNodes = Pyramid_nbN;
264       else if ( myVolumeNbNodes == 6 )
265         myFaceNbNodes = Penta_nbN;
266       else
267         myFaceNbNodes = Hexa_nbN;
268
269       // define volume orientation
270       XYZ botNormal;
271       GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z );
272       const SMDS_MeshNode* topNode = myVolumeNodes[ myVolumeNbNodes - 1 ];
273       const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ];
274       XYZ upDir (topNode->X() - botNode->X(),
275                  topNode->Y() - botNode->Y(),
276                  topNode->Z() - botNode->Z() );
277       myVolForward = ( botNormal.Dot( upDir ) < 0 );
278       break;
279     }
280     default: myVolume = 0;
281     }
282   }
283   return ( myVolume != 0 );
284 }
285
286 //=======================================================================
287 //function : GetInverseNodes
288 //purpose  : Return nodes vector of an inverse volume
289 //=======================================================================
290
291 #define SWAP_NODES(nodes,i1,i2)           \
292 {                                         \
293   const SMDS_MeshNode* tmp = nodes[ i1 ]; \
294   nodes[ i1 ] = nodes[ i2 ];              \
295   nodes[ i2 ] = tmp;                      \
296 }
297 void SMDS_VolumeTool::Inverse ()
298 {
299   if ( !myVolume ) return;
300
301   myVolForward = !myVolForward;
302   myCurFace = -1;
303
304   // inverse top and bottom faces
305   switch ( myVolumeNbNodes ) {
306   case 4:
307     SWAP_NODES( myVolumeNodes, 1, 2 );
308     break;
309   case 5:
310     SWAP_NODES( myVolumeNodes, 1, 3 );
311     break;
312   case 6:
313     SWAP_NODES( myVolumeNodes, 1, 2 );
314     SWAP_NODES( myVolumeNodes, 4, 5 );
315     break;
316   case 8:
317     SWAP_NODES( myVolumeNodes, 1, 3 );
318     SWAP_NODES( myVolumeNodes, 5, 7 );
319     break;
320   default:;
321   }
322 }
323
324 //=======================================================================
325 //function : GetSize
326 //purpose  : Return element volume
327 //=======================================================================
328
329 double SMDS_VolumeTool::GetSize() const
330 {
331   return 0;
332 }
333
334 //=======================================================================
335 //function : GetBaryCenter
336 //purpose  : 
337 //=======================================================================
338
339 bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
340 {
341   X = Y = Z = 0.;
342   if ( !myVolume )
343     return false;
344
345   for ( int i = 0; i < myVolumeNbNodes; i++ ) {
346     X += myVolumeNodes[ i ]->X();
347     Y += myVolumeNodes[ i ]->Y();
348     Z += myVolumeNodes[ i ]->Z();
349   }
350   X /= myVolumeNbNodes;
351   Y /= myVolumeNbNodes;
352   Z /= myVolumeNbNodes;
353
354   return true;
355 }
356
357 //=======================================================================
358 //function : SetExternalNormal
359 //purpose  : Node order will be so that faces normals are external
360 //=======================================================================
361
362 void SMDS_VolumeTool::SetExternalNormal ()
363 {
364   myExternalFaces = true;
365   myCurFace = -1;
366 }
367
368 //=======================================================================
369 //function : NbFaceNodes
370 //purpose  : Return number of nodes in the array of face nodes
371 //=======================================================================
372
373 int SMDS_VolumeTool::NbFaceNodes( int faceIndex )
374 {
375   if ( !setFace( faceIndex ))
376     return 0;
377   return myFaceNbNodes[ faceIndex ];
378 }
379
380 //=======================================================================
381 //function : GetFaceNodes
382 //purpose  : Return pointer to the array of face nodes.
383 //           To comfort link iteration, the array
384 //           length == NbFaceNodes( faceIndex ) + 1 and
385 //           the last node == the first one.
386 //=======================================================================
387
388 const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex )
389 {
390   if ( !setFace( faceIndex ))
391     return 0;
392   return myFaceNodes;
393 }
394
395 //=======================================================================
396 //function : GetFaceNodesIndices
397 //purpose  : Return pointer to the array of face nodes indices
398 //           To comfort link iteration, the array
399 //           length == NbFaceNodes( faceIndex ) + 1 and
400 //           the last node index == the first one.
401 //=======================================================================
402
403 const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex )
404 {
405   if ( !setFace( faceIndex ))
406     return 0;
407   return myFaceNodeIndices;
408 }
409
410 //=======================================================================
411 //function : GetFaceNodes
412 //purpose  : Return a set of face nodes.
413 //=======================================================================
414
415 bool SMDS_VolumeTool::GetFaceNodes (int                        faceIndex,
416                                     set<const SMDS_MeshNode*>& theFaceNodes )
417 {
418   if ( !setFace( faceIndex ))
419     return false;
420
421   theFaceNodes.clear();
422   int iNode, nbNode = myFaceNbNodes[ faceIndex ];
423   for ( iNode = 0; iNode < nbNode; iNode++ )
424     theFaceNodes.insert( myFaceNodes[ iNode ]);
425   
426   return true;
427 }
428
429 //=======================================================================
430 //function : IsFaceExternal
431 //purpose  : Check normal orientation of a returned face
432 //=======================================================================
433
434 bool SMDS_VolumeTool::IsFaceExternal( int faceIndex )
435 {
436   if ( myExternalFaces || !myVolume )
437     return true;
438
439   switch ( myVolumeNbNodes ) {
440   case 4:
441   case 5:
442     // only the bottom of a reversed tetrahedron can be internal
443     return ( myVolForward || faceIndex != 0 );
444   case 6:
445     // in a forward pentahedron, the top is internal, in a reversed one - bottom
446     return ( myVolForward ? faceIndex != 1 : faceIndex != 0 );
447   case 8: {
448     // in a forward hexahedron, even face normal is external, odd - internal
449     bool odd = faceIndex % 2;
450     return ( myVolForward ? !odd : odd );
451   }
452   default:;
453   }
454   return false;
455 }
456
457 //=======================================================================
458 //function : GetFaceNormal
459 //purpose  : Return a normal to a face
460 //=======================================================================
461
462 bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z)
463 {
464   if ( !setFace( faceIndex ))
465     return false;
466
467   XYZ p1 ( myFaceNodes[0] );
468   XYZ p2 ( myFaceNodes[1] );
469   XYZ p3 ( myFaceNodes[2] );
470   XYZ aVec12( p2 - p1 );
471   XYZ aVec13( p3 - p1 );
472   XYZ cross = aVec12.Crossed( aVec13 );
473
474   if ( myFaceNbNodes[ faceIndex ] == 4 ) {
475     XYZ p4 ( myFaceNodes[3] );
476     XYZ aVec14( p4 - p1 );
477     XYZ cross2 = aVec13.Crossed( aVec14 );
478     cross.x += cross2.x;
479     cross.y += cross2.y;
480     cross.z += cross2.z;    
481   }
482
483   double size = cross.Magnitude();
484   if ( size <= DBL_MIN )
485     return false;
486
487   X = cross.x / size;
488   Y = cross.y / size;
489   Z = cross.z / size;
490
491   return true;
492 }
493
494
495 //=======================================================================
496 //function : GetFaceArea
497 //purpose  : Return face area
498 //=======================================================================
499
500 double SMDS_VolumeTool::GetFaceArea( int faceIndex )
501 {
502   if ( !setFace( faceIndex ))
503     return 0;
504
505   XYZ p1 ( myFaceNodes[0] );
506   XYZ p2 ( myFaceNodes[1] );
507   XYZ p3 ( myFaceNodes[2] );
508   XYZ aVec12( p2 - p1 );
509   XYZ aVec13( p3 - p1 );
510   double area = aVec12.Crossed( aVec13 ).Magnitude() * 0.5;
511
512   if ( myFaceNbNodes[ faceIndex ] == 4 ) {
513     XYZ p4 ( myFaceNodes[3] );
514     XYZ aVec14( p4 - p1 );
515     area += aVec14.Crossed( aVec13 ).Magnitude() * 0.5;
516   }
517   return area;
518 }
519
520 //=======================================================================
521 //function : GetOppFaceIndex
522 //purpose  : Return index of the opposite face if it exists, else -1.
523 //=======================================================================
524
525 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
526 {
527   int ind = -1;
528   if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
529     switch ( myVolumeNbNodes ) {
530     case 6:
531       if ( faceIndex == 0 || faceIndex == 1 )
532         ind = 1 - faceIndex;
533       break;
534     case 8:
535       ind = faceIndex + ( faceIndex % 2 ? -1 : 1 );
536       break;
537     default:;
538     }
539   }
540   return ind;
541 }
542
543 //=======================================================================
544 //function : IsLinked
545 //purpose  : return true if theNode1 is linked with theNode2
546 //=======================================================================
547
548 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
549                                 const SMDS_MeshNode* theNode2) const
550 {
551   if ( !myVolume )
552     return false;
553
554   // find nodes indices
555   int i1 = -1, i2 = -1;
556   for ( int i = 0; i < myVolumeNbNodes; i++ ) {
557     if ( myVolumeNodes[ i ] == theNode1 )
558       i1 = i;
559     else if ( myVolumeNodes[ i ] == theNode2 )
560       i2 = i;
561   }
562   return IsLinked( i1, i2 );
563 }
564
565 //=======================================================================
566 //function : IsLinked
567 //purpose  : return true if the node with theNode1Index is linked
568 //           with the node with theNode2Index
569 //=======================================================================
570
571 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
572                                 const int theNode2Index) const
573 {
574   int minInd = theNode1Index < theNode2Index ? theNode1Index : theNode2Index;
575   int maxInd = theNode1Index < theNode2Index ? theNode2Index : theNode1Index;
576
577   if ( minInd < 0 || maxInd > myVolumeNbNodes - 1 || maxInd == minInd )
578     return false;
579
580   switch ( myVolumeNbNodes ) {
581   case 4:
582     return true;
583   case 5:
584     if ( maxInd == 4 )
585       return true;
586     switch ( maxInd - minInd ) {
587     case 1:
588     case 3: return true;
589     default:;
590     }
591     break;
592   case 6:
593     switch ( maxInd - minInd ) {
594     case 1: return minInd != 2;
595     case 2: return minInd == 0 || minInd == 3;
596     case 3: return true;
597     default:;
598     }
599     break;
600   case 8:
601     switch ( maxInd - minInd ) {
602     case 1: return minInd != 3;
603     case 3: return minInd == 0 || minInd == 4;
604     case 4: return true;
605     default:;
606     }
607     break;
608   default:;
609   }
610   return false;
611 }
612
613 //=======================================================================
614 //function : GetNodeIndex
615 //purpose  : Return an index of theNode
616 //=======================================================================
617
618 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
619 {
620   if ( myVolume ) {
621     for ( int i = 0; i < myVolumeNbNodes; i++ ) {
622       if ( myVolumeNodes[ i ] == theNode )
623         return i;
624     }
625   }
626   return -1;
627 }
628
629
630 //=======================================================================
631 //function : IsFreeFace
632 //purpose  : check that only one volume is build on the face nodes
633 //=======================================================================
634
635 bool SMDS_VolumeTool::IsFreeFace( int faceIndex )
636 {
637   const int free = true;
638   if ( !setFace( faceIndex ))
639     return !free;
640
641   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
642   int nbFaceNodes = NbFaceNodes( faceIndex );
643
644   // evaluate nb of face nodes shared by other volume
645   int maxNbShared = -1;
646   typedef map< const SMDS_MeshElement*, int > TElemIntMap;
647   TElemIntMap volNbShared;
648   TElemIntMap::iterator vNbIt;
649   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
650   {
651     const SMDS_MeshNode* n = nodes[ iNode ];
652     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator();
653     while ( eIt->more() ) {
654       const SMDS_MeshElement* elem = eIt->next();
655       if ( elem != myVolume && elem->GetType() == SMDSAbs_Volume ) {
656         int nbShared = 1;
657         vNbIt = volNbShared.find( elem );
658         if ( vNbIt == volNbShared.end() )
659           volNbShared.insert ( TElemIntMap::value_type( elem, nbShared ));
660         else
661           nbShared = ++(*vNbIt).second;
662         if ( nbShared > maxNbShared )
663           maxNbShared = nbShared;
664       }
665     }
666   }
667   if ( maxNbShared < 3 )
668     return free; // is free
669
670   // find volumes laying on the opposite side of the face
671   // and sharing all nodes
672   XYZ intNormal; // internal normal
673   GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
674   if ( IsFaceExternal( faceIndex ))
675     intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
676   XYZ p0 ( nodes[0] ), baryCenter;
677   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ )
678   {
679     int nbShared = (*vNbIt).second;
680     if ( nbShared >= 3 ) {
681       SMDS_VolumeTool volume( (*vNbIt).first );
682       volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
683       XYZ intNormal2( baryCenter - p0 );
684       if ( intNormal.Dot( intNormal2 ) < 0 )
685         continue; // opposite side
686     }
687     // remove a volume from volNbShared map
688     volNbShared.erase( vNbIt );
689   }
690   // here volNbShared contains only volumes laying on the
691   // opposite side of the face
692   if ( volNbShared.empty() )
693     return free; // is free
694
695   // check if the whole area of a face is shared
696   bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
697   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ )
698   {
699     SMDS_VolumeTool volume( (*vNbIt).first );
700     bool prevLinkShared = false;
701     int nbSharedLinks = 0;
702     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
703     {
704       bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
705       if ( linkShared )
706         nbSharedLinks++;
707       if ( linkShared && prevLinkShared &&
708           volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
709         isShared[ iNode ] = true;
710       prevLinkShared = linkShared;
711     }
712     if ( nbSharedLinks == nbFaceNodes )
713       return !free; // is not free
714     if ( nbFaceNodes == 4 ) {
715       // check traingle parts 1 & 3
716       if ( isShared[1] && isShared[3] )
717         return !free; // is not free
718       // check traingle parts 0 & 2;
719       // 0 part could not be checked in the loop; check it here
720       if ( isShared[2] && prevLinkShared &&
721           volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
722           volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
723         return !free; // is not free
724     }
725   }
726   return free;
727 }
728
729 //=======================================================================
730 //function : GetFaceIndex
731 //purpose  : Return index of a face formed by theFaceNodes
732 //=======================================================================
733
734 int SMDS_VolumeTool::GetFaceIndex( const set<const SMDS_MeshNode*>& theFaceNodes )
735 {
736   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
737     const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
738     int nbFaceNodes = NbFaceNodes( iFace );
739     set<const SMDS_MeshNode*> nodeSet;
740     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
741       nodeSet.insert( nodes[ iNode ] );
742     if ( theFaceNodes == nodeSet )
743       return iFace;
744   }
745   return -1;
746 }
747
748 //=======================================================================
749 //function : GetFaceIndex
750 //purpose  : Return index of a face formed by theFaceNodes
751 //=======================================================================
752
753 int SMDS_VolumeTool::GetFaceIndex( const set<int>& theFaceNodesIndices )
754 {
755   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
756     const int* nodes = GetFaceNodesIndices( iFace );
757     int nbFaceNodes = NbFaceNodes( iFace );
758     set<int> nodeSet;
759     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
760       nodeSet.insert( nodes[ iNode ] );
761     if ( theFaceNodesIndices == nodeSet )
762       return iFace;
763   }
764   return -1;
765 }
766
767 //=======================================================================
768 //function : setFace
769 //purpose  : 
770 //=======================================================================
771
772 bool SMDS_VolumeTool::setFace( int faceIndex )
773 {
774   if ( !myVolume )
775     return false;
776
777   if ( myCurFace == faceIndex )
778     return true;
779
780   if ( faceIndex < 0 || faceIndex >= NbFaces() )
781     return false;
782
783   // choose face node indices
784   switch ( myVolumeNbNodes ) {
785   case 4:
786     if ( myExternalFaces )
787       myFaceNodeIndices = myVolForward ? Tetra_F[ faceIndex ] : Tetra_RE[ faceIndex ];
788     else
789       myFaceNodeIndices = myVolForward ? Tetra_F[ faceIndex ] : Tetra_R[ faceIndex ];
790     break;
791   case 5:
792     if ( myExternalFaces )
793       myFaceNodeIndices = myVolForward ? Pyramid_F[ faceIndex ] : Pyramid_RE[ faceIndex ];
794     else
795       myFaceNodeIndices = myVolForward ? Pyramid_F[ faceIndex ] : Pyramid_R[ faceIndex ];
796     break;
797   case 6:
798     if ( myExternalFaces )
799       myFaceNodeIndices = myVolForward ? Penta_FE[ faceIndex ] : Penta_RE[ faceIndex ];
800     else
801       myFaceNodeIndices = myVolForward ? Penta_F[ faceIndex ] : Penta_R[ faceIndex ];
802     break;
803   case 8:
804     if ( myExternalFaces )
805       myFaceNodeIndices = myVolForward ? Hexa_FE[ faceIndex ] : Hexa_RE[ faceIndex ];
806     else
807       myFaceNodeIndices = Hexa_F[ faceIndex ];
808     break;
809   default: return false;
810   }
811
812   // set face nodes
813   int iNode, nbNode = myFaceNbNodes[ faceIndex ];
814   for ( iNode = 0; iNode < nbNode; iNode++ )
815     myFaceNodes[ iNode ] = myVolumeNodes[ myFaceNodeIndices[ iNode ]];
816   myFaceNodes[ iNode ] = myFaceNodes[ 0 ];
817
818   myCurFace = faceIndex;
819
820   return true;
821 }