Salome HOME
Update French translation file
[modules/smesh.git] / src / SMDS / SMDS_VolumeTool.cxx
index 9f35e015818776d4b5f64f2383a56492ebdd2086..745ed6454018bfd825e4970c29ca42fb4c06e2bc 100644 (file)
@@ -1,27 +1,29 @@
-// Copyright (C) 2005  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// Copyright (C) 2007-2011  CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
-// 
+//
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either 
+// License as published by the Free Software Foundation; either
 // version 2.1 of the License.
-// 
-// This library is distributed in the hope that it will be useful 
-// but WITHOUT ANY WARRANTY; without even the implied warranty of 
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Lesser General Public License for more details.
 //
-// You should have received a copy of the GNU Lesser General Public  
-// License along with this library; if not, write to the Free Software 
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
-// See http://www.salome-platform.org/
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
+
 // File      : SMDS_VolumeTool.cxx
 // Created   : Tue Jul 13 12:22:13 2004
 // Author    : Edward AGAPOV (eap)
-// Copyright : Open CASCADE
-
+//
 #ifdef _MSC_VER
 #pragma warning(disable:4786)
 #endif
 
 #include "SMDS_MeshElement.hxx"
 #include "SMDS_MeshNode.hxx"
-#include "SMDS_PolyhedralVolumeOfNodes.hxx"
+#include "SMDS_VtkVolume.hxx"
+#include "SMDS_Mesh.hxx"
 
 #include "utilities.h"
 
 #include <map>
-#include <float.h>
-#include <math.h>
+#include <limits>
+#include <cmath>
 
 using namespace std;
 
@@ -366,6 +369,8 @@ static int QuadHexa_nbN [] = { 8, 8, 8, 8, 8, 8 };
 // ========================================================
 // to perform some calculations without linkage to CASCADE
 // ========================================================
+namespace
+{
 struct XYZ {
   double x;
   double y;
@@ -374,25 +379,30 @@ struct XYZ {
   XYZ( double X, double Y, double Z ) { x = X; y = Y; z = Z; }
   XYZ( const XYZ& other )             { x = other.x; y = other.y; z = other.z; }
   XYZ( const SMDS_MeshNode* n )       { x = n->X(); y = n->Y(); z = n->Z(); }
-  XYZ operator-( const XYZ& other );
-  XYZ Crossed( const XYZ& other );
-  double Dot( const XYZ& other );
-  double Magnitude();
+  inline XYZ operator-( const XYZ& other );
+  inline XYZ Crossed( const XYZ& other );
+  inline double Dot( const XYZ& other );
+  inline double Magnitude();
+  inline double SquareMagnitude();
 };
-XYZ XYZ::operator-( const XYZ& Right ) {
+inline XYZ XYZ::operator-( const XYZ& Right ) {
   return XYZ(x - Right.x, y - Right.y, z - Right.z);
 }
-XYZ XYZ::Crossed( const XYZ& Right ) {
+inline XYZ XYZ::Crossed( const XYZ& Right ) {
   return XYZ (y * Right.z - z * Right.y,
               z * Right.x - x * Right.z,
               x * Right.y - y * Right.x);
 }
-double XYZ::Dot( const XYZ& Other ) {
+inline double XYZ::Dot( const XYZ& Other ) {
   return(x * Other.x + y * Other.y + z * Other.z);
 }
-double XYZ::Magnitude() {
+inline double XYZ::Magnitude() {
   return sqrt (x * x + y * y + z * z);
 }
+inline double XYZ::SquareMagnitude() {
+  return (x * x + y * y + z * z);
+}
+}
 
 //=======================================================================
 //function : SMDS_VolumeTool
@@ -412,6 +422,7 @@ SMDS_VolumeTool::SMDS_VolumeTool ()
        myFaceNodeIndices( NULL ),
        myFaceNodes( NULL )
 {
+  //MESSAGE("******************************************************** SMDS_VolumeToo");
 }
 
 //=======================================================================
@@ -432,6 +443,7 @@ SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume)
        myFaceNodeIndices( NULL ),
        myFaceNodes( NULL )
 {
+  //MESSAGE("******************************************************** SMDS_VolumeToo");
   Set( theVolume );
 }
 
@@ -442,14 +454,11 @@ SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume)
 
 SMDS_VolumeTool::~SMDS_VolumeTool()
 {
-  if (myVolumeNodes != NULL) {
-    delete [] myVolumeNodes;
-    myVolumeNodes = NULL;
-  }
-  if (myFaceNodes != NULL) {
-    delete [] myFaceNodes;
-    myFaceNodes = NULL;
-  }
+  if ( myVolumeNodes != NULL ) delete [] myVolumeNodes;
+  if ( myFaceNodes != NULL   ) delete [] myFaceNodes;
+
+  myFaceNodeIndices = NULL;
+  myVolumeNodes = myFaceNodes = NULL;
 }
 
 //=======================================================================
@@ -496,7 +505,7 @@ bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume)
     }
 
     if (myVolume->IsPoly()) {
-      myPolyedre = static_cast<const SMDS_PolyhedralVolumeOfNodes*>( myVolume );
+      myPolyedre = dynamic_cast<const SMDS_VtkVolume*>( myVolume );
       if (!myPolyedre) {
         MESSAGE("Warning: bad volumic element");
         return false;
@@ -515,8 +524,10 @@ bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume)
         // define volume orientation
         XYZ botNormal;
         GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z );
-        const SMDS_MeshNode* topNode = myVolumeNodes[ myVolumeNbNodes - 1 ];
         const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ];
+        int topNodeIndex = myVolume->NbCornerNodes() - 1;
+        while ( !IsLinked( 0, topNodeIndex, /*ignoreMediumNodes=*/true )) --topNodeIndex;
+        const SMDS_MeshNode* topNode = myVolumeNodes[ topNodeIndex ];
         XYZ upDir (topNode->X() - botNode->X(),
                    topNode->Y() - botNode->Y(),
                    topNode->Z() - botNode->Z() );
@@ -689,12 +700,13 @@ double SMDS_VolumeTool::GetSize() const
     if ( !myPolyedre )
       return 0.;
 
+    SMDS_Mesh *mesh = SMDS_Mesh::_meshList[myPolyedre->getMeshId()];
     // split a polyhedron into tetrahedrons
 
     SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
     XYZ baryCenter;
     me->GetBaryCenter(baryCenter.x, baryCenter.y, baryCenter.z);
-    SMDS_MeshNode bcNode ( baryCenter.x, baryCenter.y, baryCenter.z );
+    SMDS_MeshNode *bcNode = mesh->AddNode( baryCenter.x, baryCenter.y, baryCenter.z );
 
     for ( int f = 0; f < NbFaces(); ++f )
     {
@@ -704,11 +716,12 @@ double SMDS_VolumeTool::GetSize() const
         double Vn = getTetraVolume( myFaceNodes[ 0 ],
                                     myFaceNodes[ n-1 ],
                                     myFaceNodes[ n ],
-                                    bcNode );
+                                    bcNode );
 ///         cout <<"++++   " << Vn << "   nodes " <<myFaceNodes[ 0 ]->GetID() << " " <<myFaceNodes[ n-1 ]->GetID() << " " <<myFaceNodes[ n ]->GetID() << "        < " << V << endl;
         V += externalFace ? -Vn : Vn;
       }
     }
+    mesh->RemoveNode(bcNode);
   }
   else 
   {
@@ -835,6 +848,32 @@ bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
   return true;
 }
 
+//================================================================================
+/*!
+ * \brief Classify a point
+ *  \param tol - thickness of faces
+ */
+//================================================================================
+
+bool SMDS_VolumeTool::IsOut(double X, double Y, double Z, double tol) const
+{
+  // LIMITATION: for convex volumes only
+  XYZ p( X,Y,Z );
+  for ( int iF = 0; iF < myNbFaces; ++iF )
+  {
+    XYZ faceNormal;
+    if ( !GetFaceNormal( iF, faceNormal.x, faceNormal.y, faceNormal.z ))
+      continue;
+    if ( !IsFaceExternal( iF ))
+      faceNormal = XYZ() - faceNormal; // reverse
+
+    XYZ face2p( p - XYZ( myFaceNodes[0] ));
+    if ( face2p.Dot( faceNormal ) > tol )
+      return true;
+  }
+  return false;
+}
+
 //=======================================================================
 //function : SetExternalNormal
 //purpose  : Node order will be so that faces normals are external
@@ -851,7 +890,7 @@ void SMDS_VolumeTool::SetExternalNormal ()
 //purpose  : Return number of nodes in the array of face nodes
 //=======================================================================
 
-int SMDS_VolumeTool::NbFaceNodes( int faceIndex )
+int SMDS_VolumeTool::NbFaceNodes( int faceIndex ) const
 {
     if ( !setFace( faceIndex ))
       return 0;
@@ -866,7 +905,7 @@ int SMDS_VolumeTool::NbFaceNodes( int faceIndex )
 //           the last node == the first one.
 //=======================================================================
 
-const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex )
+const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex ) const
 {
   if ( !setFace( faceIndex ))
     return 0;
@@ -881,14 +920,19 @@ const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex )
 //           the last node index == the first one.
 //=======================================================================
 
-const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex )
+const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex ) const
 {
-  if (myVolume->IsPoly()) {
-    MESSAGE("Warning: attempt to obtain FaceNodesIndices of polyhedral volume");
-    return NULL;
-  }
   if ( !setFace( faceIndex ))
     return 0;
+
+  if (myVolume->IsPoly())
+  {
+    SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
+    me->myPolyIndices.resize( myFaceNbNodes + 1 );
+    me->myFaceNodeIndices = & me->myPolyIndices[0];
+    for ( int i = 0; i <= myFaceNbNodes; ++i )
+      me->myFaceNodeIndices[i] = myVolume->GetNodeIndex( myFaceNodes[i] );
+  }
   return myFaceNodeIndices;
 }
 
@@ -898,7 +942,7 @@ const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex )
 //=======================================================================
 
 bool SMDS_VolumeTool::GetFaceNodes (int                        faceIndex,
-                                    set<const SMDS_MeshNode*>& theFaceNodes )
+                                    set<const SMDS_MeshNode*>& theFaceNodes ) const
 {
   if ( !setFace( faceIndex ))
     return false;
@@ -913,10 +957,10 @@ bool SMDS_VolumeTool::GetFaceNodes (int                        faceIndex,
 
 //=======================================================================
 //function : IsFaceExternal
-//purpose  : Check normal orientation of a returned face
+//purpose  : Check normal orientation of a given face
 //=======================================================================
 
-bool SMDS_VolumeTool::IsFaceExternal( int faceIndex )
+bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) const
 {
   if ( myExternalFaces || !myVolume )
     return true;
@@ -958,7 +1002,7 @@ bool SMDS_VolumeTool::IsFaceExternal( int faceIndex )
 //purpose  : Return a normal to a face
 //=======================================================================
 
-bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z)
+bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z) const
 {
   if ( !setFace( faceIndex ))
     return false;
@@ -981,7 +1025,7 @@ bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, doub
   }
 
   double size = cross.Magnitude();
-  if ( size <= DBL_MIN )
+  if ( size <= numeric_limits<double>::min() )
     return false;
 
   X = cross.x / size;
@@ -991,12 +1035,33 @@ bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, doub
   return true;
 }
 
+//================================================================================
+/*!
+ * \brief Return barycenter of a face
+ */
+//================================================================================
+
+bool SMDS_VolumeTool::GetFaceBaryCenter (int faceIndex, double & X, double & Y, double & Z) const
+{
+  if ( !setFace( faceIndex ))
+    return false;
+
+  X = Y = Z = 0.0;
+  for ( int i = 0; i < myFaceNbNodes; ++i )
+  {
+    X += myFaceNodes[i]->X() / myFaceNbNodes;
+    Y += myFaceNodes[i]->Y() / myFaceNbNodes;
+    Z += myFaceNodes[i]->Z() / myFaceNbNodes;
+  }
+  return true;
+}
+
 //=======================================================================
 //function : GetFaceArea
 //purpose  : Return face area
 //=======================================================================
 
-double SMDS_VolumeTool::GetFaceArea( int faceIndex )
+double SMDS_VolumeTool::GetFaceArea( int faceIndex ) const
 {
   if (myVolume->IsPoly()) {
     MESSAGE("Warning: attempt to obtain area of a face of polyhedral volume");
@@ -1052,10 +1117,12 @@ int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
 //=======================================================================
 //function : IsLinked
 //purpose  : return true if theNode1 is linked with theNode2
+// If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
 //=======================================================================
 
 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
-                                const SMDS_MeshNode* theNode2) const
+                                const SMDS_MeshNode* theNode2,
+                                const bool           theIgnoreMediumNodes) const
 {
   if ( !myVolume )
     return false;
@@ -1102,25 +1169,50 @@ bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
 //function : IsLinked
 //purpose  : return true if the node with theNode1Index is linked
 //           with the node with theNode2Index
+// If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
 //=======================================================================
 
 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
-                                const int theNode2Index) const
+                                const int theNode2Index,
+                                bool      theIgnoreMediumNodes) const
 {
   if ( myVolume->IsPoly() ) {
     return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
   }
 
-  int minInd = theNode1Index < theNode2Index ? theNode1Index : theNode2Index;
-  int maxInd = theNode1Index < theNode2Index ? theNode2Index : theNode1Index;
+  int minInd = min( theNode1Index, theNode2Index );
+  int maxInd = max( theNode1Index, theNode2Index );
 
   if ( minInd < 0 || maxInd > myVolumeNbNodes - 1 || maxInd == minInd )
     return false;
 
-  switch ( myVolumeNbNodes ) {
-  case 4:
+  SMDSAbs_EntityType type = myVolume->GetEntityType();
+  if ( myVolume->IsQuadratic() )
+  {
+    int firstMediumInd = myVolume->NbCornerNodes();
+    if ( minInd >= firstMediumInd )
+      return false; // medium nodes are not linked
+    if ( maxInd < firstMediumInd ) // both nodes are corners
+    {
+      if ( theIgnoreMediumNodes )
+        type = SMDSAbs_EntityType( int(type)-1 ); // check linkage of corner nodes
+      else
+        return false; // corner nodes are not linked directly in a quadratic cell
+    }
+  }
+
+  switch ( type ) {
+  case SMDSEntity_Tetra:
     return true;
-  case 5:
+  case SMDSEntity_Hexa:
+    switch ( maxInd - minInd ) {
+    case 1: return minInd != 3;
+    case 3: return minInd == 0 || minInd == 4;
+    case 4: return true;
+    default:;
+    }
+    break;
+  case SMDSEntity_Pyramid:
     if ( maxInd == 4 )
       return true;
     switch ( maxInd - minInd ) {
@@ -1129,7 +1221,7 @@ bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
     default:;
     }
     break;
-  case 6:
+  case SMDSEntity_Penta:
     switch ( maxInd - minInd ) {
     case 1: return minInd != 2;
     case 2: return minInd == 0 || minInd == 3;
@@ -1137,15 +1229,7 @@ bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
     default:;
     }
     break;
-  case 8:
-    switch ( maxInd - minInd ) {
-    case 1: return minInd != 3;
-    case 3: return minInd == 0 || minInd == 4;
-    case 4: return true;
-    default:;
-    }
-    break;
-  case 10:
+  case SMDSEntity_Quad_Tetra:
     {
       switch ( minInd ) {
       case 0: if( maxInd==4 ||  maxInd==6 ||  maxInd==7 ) return true;
@@ -1156,7 +1240,22 @@ bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
       }
       break;
     }
-  case 13:
+  case SMDSEntity_Quad_Hexa:
+    {
+      switch ( minInd ) {
+      case 0: if( maxInd==8 ||  maxInd==11 ||  maxInd==16 ) return true;
+      case 1: if( maxInd==8 ||  maxInd==9 ||  maxInd==17 ) return true;
+      case 2: if( maxInd==9 ||  maxInd==10 ||  maxInd==18 ) return true;
+      case 3: if( maxInd==10 ||  maxInd==11 ||  maxInd==19 ) return true;
+      case 4: if( maxInd==12 ||  maxInd==15 ||  maxInd==16 ) return true;
+      case 5: if( maxInd==12 ||  maxInd==13 ||  maxInd==17 ) return true;
+      case 6: if( maxInd==13 ||  maxInd==14 ||  maxInd==18 ) return true;
+      case 7: if( maxInd==14 ||  maxInd==15 ||  maxInd==19 ) return true;
+      default:;
+      }
+      break;
+    }
+  case SMDSEntity_Quad_Pyramid:
     {
       switch ( minInd ) {
       case 0: if( maxInd==5 ||  maxInd==8 ||  maxInd==9 ) return true;
@@ -1168,7 +1267,7 @@ bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
       }
       break;
     }
-  case 15:
+  case SMDSEntity_Quad_Penta:
     {
       switch ( minInd ) {
       case 0: if( maxInd==6 ||  maxInd==8 ||  maxInd==12 ) return true;
@@ -1181,21 +1280,6 @@ bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
       }
       break;
     }
-  case 20:
-    {
-      switch ( minInd ) {
-      case 0: if( maxInd==8 ||  maxInd==11 ||  maxInd==16 ) return true;
-      case 1: if( maxInd==8 ||  maxInd==9 ||  maxInd==17 ) return true;
-      case 2: if( maxInd==9 ||  maxInd==10 ||  maxInd==18 ) return true;
-      case 3: if( maxInd==10 ||  maxInd==11 ||  maxInd==19 ) return true;
-      case 4: if( maxInd==12 ||  maxInd==15 ||  maxInd==16 ) return true;
-      case 5: if( maxInd==12 ||  maxInd==13 ||  maxInd==17 ) return true;
-      case 6: if( maxInd==13 ||  maxInd==14 ||  maxInd==18 ) return true;
-      case 7: if( maxInd==14 ||  maxInd==15 ||  maxInd==19 ) return true;
-      default:;
-      }
-      break;
-    }
   default:;
   }
   return false;
@@ -1217,20 +1301,121 @@ int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
   return -1;
 }
 
-//=======================================================================
-//function : IsFreeFace
-//purpose  : check that only one volume is build on the face nodes
-//=======================================================================
+//================================================================================
+/*!
+ * \brief Fill vector with boundary faces existing in the mesh
+  * \param faces - vector of found nodes
+  * \retval int - nb of found faces
+ */
+//================================================================================
+
+int SMDS_VolumeTool::GetAllExistingFaces(vector<const SMDS_MeshElement*> & faces) const
+{
+  faces.clear();
+  faces.reserve( NbFaces() );
+  for ( int iF = 0; iF < NbFaces(); ++iF ) {
+    const SMDS_MeshFace* face = 0;
+    const SMDS_MeshNode** nodes = GetFaceNodes( iF );
+    switch ( NbFaceNodes( iF )) {
+    case 3:
+      face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2] ); break;
+    case 4:
+      face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
+    case 6:
+      face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2],
+                                  nodes[3], nodes[4], nodes[5]); break;
+    case 8:
+      face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3],
+                                  nodes[4], nodes[5], nodes[6], nodes[7]); break;
+    }
+    if ( face )
+      faces.push_back( face );
+  }
+  return faces.size();
+}
+
+
+//================================================================================
+/*!
+ * \brief Fill vector with boundary edges existing in the mesh
+  * \param edges - vector of found edges
+  * \retval int - nb of found faces
+ */
+//================================================================================
+
+int SMDS_VolumeTool::GetAllExistingEdges(vector<const SMDS_MeshElement*> & edges) const
+{
+  edges.clear();
+  edges.reserve( myVolumeNbNodes * 2 );
+  for ( int i = 0; i < myVolumeNbNodes; ++i ) {
+    for ( int j = i + 1; j < myVolumeNbNodes; ++j ) {
+      if ( IsLinked( i, j )) {
+        const SMDS_MeshElement* edge =
+          SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] );
+        if ( edge )
+          edges.push_back( edge );
+      }
+    }
+  }
+  return edges.size();
+}
+
+//================================================================================
+/*!
+ * \brief Return minimal square distance between connected corner nodes
+ */
+//================================================================================
+
+double SMDS_VolumeTool::MinLinearSize2() const
+{
+  double minSize = 1e+100;
+  int iQ = myVolume->IsQuadratic() ? 2 : 1;
+
+  // store current face data
+  int curFace = myCurFace, nbN = myFaceNbNodes;
+  int* ind = myFaceNodeIndices;
+  myFaceNodeIndices = NULL;
+  const SMDS_MeshNode** nodes = myFaceNodes;
+  myFaceNodes = NULL;
+  
+  // it seems that compute distance twice is faster than organization of a sole computing
+  myCurFace = -1;
+  for ( int iF = 0; iF < myNbFaces; ++iF )
+  {
+    setFace( iF );
+    for ( int iN = 0; iN < myFaceNbNodes; iN += iQ )
+    {
+      XYZ n1( myFaceNodes[ iN ]);
+      XYZ n2( myFaceNodes[(iN + iQ) % myFaceNbNodes]);
+      minSize = std::min( minSize, (n1 - n2).SquareMagnitude());
+    }
+  }
+  // restore current face data
+  myCurFace = curFace;
+  myFaceNbNodes = nbN;
+  myFaceNodeIndices = ind;
+  delete [] myFaceNodes; myFaceNodes = nodes;
+
+  return minSize;
+}
+
+//================================================================================
+/*!
+ * \brief check that only one volume is build on the face nodes
+ *
+ * If a face is shared by one of <ignoreVolumes>, it is considered free
+ */
+//================================================================================
 
-bool SMDS_VolumeTool::IsFreeFace( int faceIndex )
+bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
 {
-  const int free = true;
+  const bool isFree = true;
 
   if (!setFace( faceIndex ))
-    return !free;
+    return !isFree;
 
   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
-  int nbFaceNodes = myFaceNbNodes;
+  const int nbFaceNodes = myFaceNbNodes;
 
   // evaluate nb of face nodes shared by other volume
   int maxNbShared = -1;
@@ -1239,25 +1424,19 @@ bool SMDS_VolumeTool::IsFreeFace( int faceIndex )
   TElemIntMap::iterator vNbIt;
   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
     const SMDS_MeshNode* n = nodes[ iNode ];
-    SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator();
+    SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
     while ( eIt->more() ) {
       const SMDS_MeshElement* elem = eIt->next();
-      if ( elem != myVolume && elem->GetType() == SMDSAbs_Volume ) {
-        int nbShared = 1;
-        vNbIt = volNbShared.find( elem );
-        if ( vNbIt == volNbShared.end() ) {
-          volNbShared.insert ( TElemIntMap::value_type( elem, nbShared ));
-        }
-        else {
-          nbShared = ++(*vNbIt).second;
-        }
-        if ( nbShared > maxNbShared )
-          maxNbShared = nbShared;
+      if ( elem != myVolume ) {
+        vNbIt = volNbShared.insert( make_pair( elem, 0 )).first;
+        (*vNbIt).second++;
+        if ( vNbIt->second > maxNbShared )
+          maxNbShared = vNbIt->second;
       }
     }
   }
   if ( maxNbShared < 3 )
-    return free; // is free
+    return isFree; // is free
 
   // find volumes laying on the opposite side of the face
   // and sharing all nodes
@@ -1266,55 +1445,81 @@ bool SMDS_VolumeTool::IsFreeFace( int faceIndex )
   if ( IsFaceExternal( faceIndex ))
     intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
   XYZ p0 ( nodes[0] ), baryCenter;
-  for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
-    int nbShared = (*vNbIt).second;
+  for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end();  ) {
+    const int& nbShared = (*vNbIt).second;
     if ( nbShared >= 3 ) {
       SMDS_VolumeTool volume( (*vNbIt).first );
       volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
       XYZ intNormal2( baryCenter - p0 );
-      if ( intNormal.Dot( intNormal2 ) < 0 )
-        continue; // opposite side
+      if ( intNormal.Dot( intNormal2 ) < 0 ) {
+        // opposite side
+        if ( nbShared >= nbFaceNodes )
+        {
+          // a volume shares the whole facet
+          if ( otherVol ) *otherVol = vNbIt->first;
+          return !isFree; 
+        }
+        ++vNbIt;
+        continue;
+      }
     }
     // remove a volume from volNbShared map
-    volNbShared.erase( vNbIt );
+    volNbShared.erase( vNbIt++ );
   }
 
-  // here volNbShared contains only volumes laying on the
-  // opposite side of the face
-  if ( volNbShared.empty() ) {
-    return free; // is free
+  // here volNbShared contains only volumes laying on the opposite side of
+  // the face and sharing 3 or more but not all face nodes with myVolume
+  if ( volNbShared.size() < 2 ) {
+    return isFree; // is free
   }
 
   // check if the whole area of a face is shared
-  bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
-  for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
-    SMDS_VolumeTool volume( (*vNbIt).first );
-    bool prevLinkShared = false;
-    int nbSharedLinks = 0;
-    for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
-      bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
-      if ( linkShared )
-        nbSharedLinks++;
-      if ( linkShared && prevLinkShared &&
-          volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
-        isShared[ iNode ] = true;
-      prevLinkShared = linkShared;
-    }
-    if ( nbSharedLinks == nbFaceNodes )
-      return !free; // is not free
-    if ( nbFaceNodes == 4 ) {
-      // check traingle parts 1 & 3
-      if ( isShared[1] && isShared[3] )
-        return !free; // is not free
-      // check triangle parts 0 & 2;
-      // 0 part could not be checked in the loop; check it here
-      if ( isShared[2] && prevLinkShared &&
-          volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
-          volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
-        return !free; // is not free
-    }
+  for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
+  {
+    const SMDS_MeshNode* n = nodes[ iNode ];
+    // check if n is shared by one of volumes of volNbShared
+    bool isShared = false;
+    SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
+    while ( eIt->more() && !isShared )
+      isShared = volNbShared.count( eIt->next() );
+    if ( !isShared )
+      return isFree;
   }
-  return free;
+  if ( otherVol ) *otherVol = volNbShared.begin()->first;
+  return !isFree;
+
+//   if ( !myVolume->IsPoly() )
+//   {
+//     bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
+//     for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
+//       SMDS_VolumeTool volume( (*vNbIt).first );
+//       bool prevLinkShared = false;
+//       int nbSharedLinks = 0;
+//       for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
+//         bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
+//         if ( linkShared )
+//           nbSharedLinks++;
+//         if ( linkShared && prevLinkShared &&
+//              volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
+//           isShared[ iNode ] = true;
+//         prevLinkShared = linkShared;
+//       }
+//       if ( nbSharedLinks == nbFaceNodes )
+//         return !free; // is not free
+//       if ( nbFaceNodes == 4 ) {
+//         // check traingle parts 1 & 3
+//         if ( isShared[1] && isShared[3] )
+//           return !free; // is not free
+//         // check triangle parts 0 & 2;
+//         // 0 part could not be checked in the loop; check it here
+//         if ( isShared[2] && prevLinkShared &&
+//              volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
+//              volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
+//           return !free; // is not free
+//       }
+//     }
+//   }
+//  return free;
 }
 
 //=======================================================================
@@ -1322,7 +1527,7 @@ bool SMDS_VolumeTool::IsFreeFace( int faceIndex )
 //purpose  : Return index of a face formed by theFaceNodes
 //=======================================================================
 
-int SMDS_VolumeTool::GetFaceIndex( const set<const SMDS_MeshNode*>& theFaceNodes )
+int SMDS_VolumeTool::GetFaceIndex( const set<const SMDS_MeshNode*>& theFaceNodes ) const
 {
   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
     const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
@@ -1360,7 +1565,7 @@ int SMDS_VolumeTool::GetFaceIndex( const set<const SMDS_MeshNode*>& theFaceNodes
 //purpose  : 
 //=======================================================================
 
-bool SMDS_VolumeTool::setFace( int faceIndex )
+bool SMDS_VolumeTool::setFace( int faceIndex ) const
 {
   if ( !myVolume )
     return false;
@@ -1384,24 +1589,24 @@ bool SMDS_VolumeTool::setFace( int faceIndex )
       return false;
     }
 
-    // check orientation
-    bool isGoodOri = true;
-    if (myExternalFaces)
-      isGoodOri = IsFaceExternal( faceIndex );
-
     // set face nodes
     int iNode;
     myFaceNbNodes = myPolyedre->NbFaceNodes(faceIndex + 1);
     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
-    if (isGoodOri) {
-      for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
-        myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, iNode + 1);
-    } else {
-      for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
-        myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, myFaceNbNodes - iNode);
-    }
+    for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
+      myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, iNode + 1);
     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ]; // last = first
 
+    // check orientation
+    if (myExternalFaces)
+    {
+      myCurFace = faceIndex; // avoid infinite recursion in IsFaceExternal()
+      myExternalFaces = false; // force normal computation by IsFaceExternal()
+      if ( !IsFaceExternal( faceIndex ))
+        for ( int i = 0, j = myFaceNbNodes; i < j; ++i, --j )
+          std::swap( myFaceNodes[i], myFaceNodes[j] );
+      myExternalFaces = true;
+    }
   }
   else {
     // choose face node indices
@@ -1518,6 +1723,31 @@ int SMDS_VolumeTool::NbFaces( VolumeType type )
   }
 }
 
+//================================================================================
+/*!
+ * \brief Useful to know nb of corner nodes of a quadratic volume
+  * \param type - volume type
+  * \retval int - nb of corner nodes
+ */
+//================================================================================
+
+int SMDS_VolumeTool::NbCornerNodes(VolumeType type)
+{
+  switch ( type ) {
+  case TETRA     :
+  case QUAD_TETRA: return 4;
+  case PYRAM     :
+  case QUAD_PYRAM: return 5;
+  case PENTA     :
+  case QUAD_PENTA: return 6;
+  case HEXA      :
+  case QUAD_HEXA : return 8;
+  default:    return 0;
+  }
+  return 0;
+}
+  // 
+
 //=======================================================================
 //function : GetFaceNodesIndices
 //purpose  : Return the array of face nodes indices
@@ -1566,3 +1796,22 @@ int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
   return 0;
 }
 
+//=======================================================================
+//function : Get
+//purpose  : return element
+//=======================================================================
+
+const SMDS_MeshVolume* SMDS_VolumeTool::Element() const
+{
+  return static_cast<const SMDS_MeshVolume*>( myVolume );
+}
+
+//=======================================================================
+//function : ID
+//purpose  : return element ID
+//=======================================================================
+
+int SMDS_VolumeTool::ID() const
+{
+  return myVolume ? myVolume->GetID() : 0;
+}