diff --git a/doc/salome/gui/SMESH/images/revolutionsn1.png b/doc/salome/gui/SMESH/images/revolutionsn1.png
index 6564c7fc6..3e7a7a1aa 100644
Binary files a/doc/salome/gui/SMESH/images/revolutionsn1.png and b/doc/salome/gui/SMESH/images/revolutionsn1.png differ
diff --git a/doc/salome/gui/SMESH/images/revolutionsn2.png b/doc/salome/gui/SMESH/images/revolutionsn2.png
index 67a673144..064dd6c2d 100644
Binary files a/doc/salome/gui/SMESH/images/revolutionsn2.png and b/doc/salome/gui/SMESH/images/revolutionsn2.png differ
diff --git a/doc/salome/gui/SMESH/input/extrusion.doc b/doc/salome/gui/SMESH/input/extrusion.doc
index af390a4b1..56cb32f4c 100644
--- a/doc/salome/gui/SMESH/input/extrusion.doc
+++ b/doc/salome/gui/SMESH/input/extrusion.doc
@@ -48,6 +48,11 @@ The following dialog will appear:
 
 <li>In this dialog:
 <ul>
+  <li>Use \a Selection button to specify what you are going to
+    select at a given moment, \b Nodes, \b Edges or \b Faces.
+    \image html image120.png
+    <center><em>"Selection" button</em></center>
+  </li>
   <li>Specify \b Nodes, \b Edges and \b Faces, which will be extruded, by one
     of following means:
     <ul>
diff --git a/doc/salome/gui/SMESH/input/revolution.doc b/doc/salome/gui/SMESH/input/revolution.doc
index bd11f5a1b..f0d4194df 100644
--- a/doc/salome/gui/SMESH/input/revolution.doc
+++ b/doc/salome/gui/SMESH/input/revolution.doc
@@ -68,16 +68,16 @@ The following dialog will appear:
       <ul> 
         <li> <b>Angle by Step</b> - the elements are revolved by the
           specified angle at each step (i.e. for Angle=30 and Number of
-          Steps=2, the elements will be extruded by 30 degrees twice for a
-          total of 30*2=60)
-\image html revolutionsn2.png "Example of Revolution with Angle by Step"
+          Steps=3, the elements will be extruded by 30 degrees twice for a
+          total of 30*3=90)
+\image html revolutionsn2.png "Example of Revolution with Angle by Step. Angle=30 and Number of Steps=3"
         </li>
         <li> <b>Total Angle</b> - the elements are revolved by the
           specified angle only once and the number of steps defines the
-          number of iterations (i.e. for Angle=30 and Number of Steps=2,
-          the elements will be revolved by 30/2=15 degrees twice for a
+          number of iterations (i.e. for Angle=30 and Number of Steps=3,
+          the elements will be revolved by 30/3=10 degrees twice for a
           total of 30).
-\image html revolutionsn1.png "Example of Revolution with Total Angle"
+\image html revolutionsn1.png "Example of Revolution with Total Angle. Angle=30 and Number of Steps=3"
         </li>
       </ul>
     </li>
diff --git a/src/SMDS/SMDS_Mesh.hxx b/src/SMDS/SMDS_Mesh.hxx
index bc3c74012..7c728d6f9 100644
--- a/src/SMDS/SMDS_Mesh.hxx
+++ b/src/SMDS/SMDS_Mesh.hxx
@@ -58,7 +58,7 @@
 
 #define MYASSERT(val) if (!(val)) throw SALOME_Exception(LOCALIZED("assertion not verified"));
 
-class SMDS_EXPORT SMDS_Mesh:public SMDS_MeshObject
+class SMDS_EXPORT SMDS_Mesh : public SMDS_MeshObject
 {
 public:
   friend class SMDS_MeshIDFactory;
diff --git a/src/SMDS/SMDS_VolumeTool.cxx b/src/SMDS/SMDS_VolumeTool.cxx
index b86dc16c3..8296de035 100644
--- a/src/SMDS/SMDS_VolumeTool.cxx
+++ b/src/SMDS/SMDS_VolumeTool.cxx
@@ -40,6 +40,8 @@
 #include <map>
 #include <limits>
 #include <cmath>
+#include <numeric>
+#include <algorithm>
 
 using namespace std;
 
@@ -396,7 +398,7 @@ inline double XYZ::SquareMagnitude() {
   SMDS_VolumeTool::VolumeType quadToLinear(SMDS_VolumeTool::VolumeType quadType)
   {
     SMDS_VolumeTool::VolumeType linType = SMDS_VolumeTool::VolumeType( int(quadType)-4 );
-    const int nbCornersByQuad = SMDS_VolumeTool::NbCornerNodes( quadType );
+    const int           nbCornersByQuad = SMDS_VolumeTool::NbCornerNodes( quadType );
     if ( SMDS_VolumeTool::NbCornerNodes( linType ) == nbCornersByQuad )
       return linType;
 
@@ -410,14 +412,33 @@ inline double XYZ::SquareMagnitude() {
 
 } // namespace
 
+//================================================================================
+/*!
+ * \brief Saver/restorer of a SMDS_VolumeTool::myCurFace
+ */
+//================================================================================
+
+struct SMDS_VolumeTool::SaveFacet
+{
+  SMDS_VolumeTool::Facet  mySaved;
+  SMDS_VolumeTool::Facet& myToRestore;
+  SaveFacet( SMDS_VolumeTool::Facet& facet ): myToRestore( facet )
+  {
+    mySaved = facet;
+  }
+  ~SaveFacet()
+  {
+    if ( myToRestore.myIndex != mySaved.myIndex )
+      myToRestore = mySaved;
+  }
+};
+
 //=======================================================================
 //function : SMDS_VolumeTool
-//purpose  : 
+//purpose  :
 //=======================================================================
 
 SMDS_VolumeTool::SMDS_VolumeTool ()
-  : myVolumeNodes( NULL ),
-    myFaceNodes( NULL )
 {
   Set( 0 );
 }
@@ -429,8 +450,6 @@ SMDS_VolumeTool::SMDS_VolumeTool ()
 
 SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume,
                                   const bool              ignoreCentralNodes)
-  : myVolumeNodes( NULL ),
-    myFaceNodes( NULL )
 {
   Set( theVolume, ignoreCentralNodes );
 }
@@ -442,11 +461,7 @@ SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume,
 
 SMDS_VolumeTool::~SMDS_VolumeTool()
 {
-  if ( myVolumeNodes != NULL ) delete [] myVolumeNodes;
-  if ( myFaceNodes != NULL   ) delete [] myFaceNodes;
-
-  myFaceNodeIndices = NULL;
-  myVolumeNodes = myFaceNodes = NULL;
+  myCurFace.myNodeIndices = NULL;
 }
 
 //=======================================================================
@@ -464,42 +479,37 @@ bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume,
 
   myVolForward = true;
   myNbFaces = 0;
-  myVolumeNbNodes = 0;
-  if (myVolumeNodes != NULL) {
-    delete [] myVolumeNodes;
-    myVolumeNodes = NULL;
-  }
+  myVolumeNodes.clear();
   myPolyIndices.clear();
+  myPolyQuantities.clear();
+  myPolyFacetOri.clear();
+  myFwdLinks.clear();
 
   myExternalFaces = false;
 
   myAllFacesNodeIndices_F  = 0;
-  //myAllFacesNodeIndices_FE = 0;
   myAllFacesNodeIndices_RE = 0;
   myAllFacesNbNodes        = 0;
 
-  myCurFace = -1;
-  myFaceNbNodes = 0;
-  myFaceNodeIndices = NULL;
-  if (myFaceNodes != NULL) {
-    delete [] myFaceNodes;
-    myFaceNodes = NULL;
-  }
+  myCurFace.myIndex = -1;
+  myCurFace.myNodeIndices = NULL;
+  myCurFace.myNodes.clear();
 
   // set volume data
   if ( !theVolume || theVolume->GetType() != SMDSAbs_Volume )
     return false;
 
   myVolume = theVolume;
-  if (myVolume->IsPoly())
-    myPolyedre = dynamic_cast<const SMDS_VtkVolume*>( myVolume );
-
   myNbFaces = theVolume->NbFaces();
-  myVolumeNbNodes = theVolume->NbNodes();
+  if ( myVolume->IsPoly() )
+  {
+    myPolyedre = dynamic_cast<const SMDS_VtkVolume*>( myVolume );
+    myPolyFacetOri.resize( myNbFaces, 0 );
+  }
 
   // set nodes
   int iNode = 0;
-  myVolumeNodes = new const SMDS_MeshNode* [myVolumeNbNodes];
+  myVolumeNodes.resize( myVolume->NbNodes() );
   SMDS_ElemIteratorPtr nodeIt = myVolume->nodesIterator();
   while ( nodeIt->more() )
     myVolumeNodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
@@ -512,18 +522,19 @@ bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume,
   {
     // define volume orientation
     XYZ botNormal;
-    GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z );
-    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() );
-    myVolForward = ( botNormal.Dot( upDir ) < 0 );
-
+    if ( GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z ))
+    {
+      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() );
+      myVolForward = ( botNormal.Dot( upDir ) < 0 );
+    }
     if ( !myVolForward )
-      myCurFace = -1; // previous setFace(0) didn't take myVolForward into account
+      myCurFace.myIndex = -1; // previous setFace(0) didn't take myVolForward into account
   }
   return true;
 }
@@ -549,10 +560,10 @@ void SMDS_VolumeTool::Inverse ()
   }
 
   myVolForward = !myVolForward;
-  myCurFace = -1;
+  myCurFace.myIndex = -1;
 
   // inverse top and bottom faces
-  switch ( myVolumeNbNodes ) {
+  switch ( myVolumeNodes.size() ) {
   case 4:
     SWAP_NODES( myVolumeNodes, 1, 2 );
     break;
@@ -626,7 +637,7 @@ SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetVolumeType() const
   if ( myPolyedre )
     return POLYHEDA;
 
-  switch( myVolumeNbNodes ) {
+  switch( myVolumeNodes.size() ) {
   case 4: return TETRA;
   case 5: return PYRAM;
   case 6: return PENTA;
@@ -653,28 +664,18 @@ static double getTetraVolume(const SMDS_MeshNode* n1,
                              const SMDS_MeshNode* n3,
                              const SMDS_MeshNode* n4)
 {
-  double X1 = n1->X();
-  double Y1 = n1->Y();
-  double Z1 = n1->Z();
+  double p1[3], p2[3], p3[3], p4[3];
+  n1->GetXYZ( p1 );
+  n2->GetXYZ( p2 );
+  n3->GetXYZ( p3 );
+  n4->GetXYZ( p4 );
 
-  double X2 = n2->X();
-  double Y2 = n2->Y();
-  double Z2 = n2->Z();
-
-  double X3 = n3->X();
-  double Y3 = n3->Y();
-  double Z3 = n3->Z();
-
-  double X4 = n4->X();
-  double Y4 = n4->Y();
-  double Z4 = n4->Z();
-
-  double Q1 = -(X1-X2)*(Y3*Z4-Y4*Z3);
-  double Q2 =  (X1-X3)*(Y2*Z4-Y4*Z2);
-  double R1 = -(X1-X4)*(Y2*Z3-Y3*Z2);
-  double R2 = -(X2-X3)*(Y1*Z4-Y4*Z1);
-  double S1 =  (X2-X4)*(Y1*Z3-Y3*Z1);
-  double S2 = -(X3-X4)*(Y1*Z2-Y2*Z1);
+  double Q1 = -(p1[ 0 ]-p2[ 0 ])*(p3[ 1 ]*p4[ 2 ]-p4[ 1 ]*p3[ 2 ]);
+  double Q2 =  (p1[ 0 ]-p3[ 0 ])*(p2[ 1 ]*p4[ 2 ]-p4[ 1 ]*p2[ 2 ]);
+  double R1 = -(p1[ 0 ]-p4[ 0 ])*(p2[ 1 ]*p3[ 2 ]-p3[ 1 ]*p2[ 2 ]);
+  double R2 = -(p2[ 0 ]-p3[ 0 ])*(p1[ 1 ]*p4[ 2 ]-p4[ 1 ]*p1[ 2 ]);
+  double S1 =  (p2[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p3[ 2 ]-p3[ 1 ]*p1[ 2 ]);
+  double S2 = -(p3[ 0 ]-p4[ 0 ])*(p1[ 1 ]*p2[ 2 ]-p2[ 1 ]*p1[ 2 ]);
 
   return (Q1+Q2+R1+R2+S1+S2)/6.0;
 }
@@ -697,23 +698,21 @@ double SMDS_VolumeTool::GetSize() const
 
     // split a polyhedron into tetrahedrons
 
-    int saveCurFace = myCurFace;
+    SaveFacet savedFacet( myCurFace );
     SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
     for ( int f = 0; f < NbFaces(); ++f )
     {
       me->setFace( f );
-      XYZ area (0,0,0), p1( myFaceNodes[0] );
-      for ( int n = 0; n < myFaceNbNodes; ++n )
+      XYZ area (0,0,0), p1( myCurFace.myNodes[0] );
+      for ( int n = 0; n < myCurFace.myNbNodes; ++n )
       {
-        XYZ p2( myFaceNodes[ n+1 ]);
+        XYZ p2( myCurFace.myNodes[ n+1 ]);
         area = area + p1.Crossed( p2 );
         p1 = p2;
       }
       V += p1.Dot( area );
     }
     V /= 6;
-    if ( saveCurFace > -1 && saveCurFace != myCurFace )
-      me->setFace( myCurFace );
   }
   else 
   {
@@ -844,14 +843,14 @@ bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
   if ( !myVolume )
     return false;
 
-  for ( int i = 0; i < myVolumeNbNodes; i++ ) {
+  for ( int i = 0; i < myVolumeNodes.size(); i++ ) {
     X += myVolumeNodes[ i ]->X();
     Y += myVolumeNodes[ i ]->Y();
     Z += myVolumeNodes[ i ]->Z();
   }
-  X /= myVolumeNbNodes;
-  Y /= myVolumeNbNodes;
-  Z /= myVolumeNbNodes;
+  X /= myVolumeNodes.size();
+  Y /= myVolumeNodes.size();
+  Z /= myVolumeNodes.size();
 
   return true;
 }
@@ -875,7 +874,7 @@ bool SMDS_VolumeTool::IsOut(double X, double Y, double Z, double tol) const
     if ( !IsFaceExternal( iF ))
       faceNormal = XYZ() - faceNormal; // reverse
 
-    XYZ face2p( p - XYZ( myFaceNodes[0] ));
+    XYZ face2p( p - XYZ( myCurFace.myNodes[0] ));
     if ( face2p.Dot( faceNormal ) > tol )
       return true;
   }
@@ -890,7 +889,7 @@ bool SMDS_VolumeTool::IsOut(double X, double Y, double Z, double tol) const
 void SMDS_VolumeTool::SetExternalNormal ()
 {
   myExternalFaces = true;
-  myCurFace = -1;
+  myCurFace.myIndex = -1;
 }
 
 //=======================================================================
@@ -900,9 +899,9 @@ void SMDS_VolumeTool::SetExternalNormal ()
 
 int SMDS_VolumeTool::NbFaceNodes( int faceIndex ) const
 {
-    if ( !setFace( faceIndex ))
-      return 0;
-    return myFaceNbNodes;
+  if ( !setFace( faceIndex ))
+    return 0;
+  return myCurFace.myNbNodes;
 }
 
 //=======================================================================
@@ -917,7 +916,7 @@ const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex ) const
 {
   if ( !setFace( faceIndex ))
     return 0;
-  return myFaceNodes;
+  return &myCurFace.myNodes[0];
 }
 
 //=======================================================================
@@ -933,15 +932,7 @@ const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex ) const
   if ( !setFace( faceIndex ))
     return 0;
 
-  if (myPolyedre)
-  {
-    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;
+  return myCurFace.myNodeIndices;
 }
 
 //=======================================================================
@@ -956,11 +947,44 @@ bool SMDS_VolumeTool::GetFaceNodes (int                        faceIndex,
     return false;
 
   theFaceNodes.clear();
-  theFaceNodes.insert( myFaceNodes, myFaceNodes + myFaceNbNodes );
+  theFaceNodes.insert( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
 
   return true;
 }
 
+namespace
+{
+  struct NLink : public std::pair<int,int>
+  {
+    int myOri;
+    NLink(const SMDS_MeshNode* n1=0, const SMDS_MeshNode* n2=0, int ori=1 )
+    {
+      if ( n1 )
+      {
+        if (( myOri = ( n1->GetID() < n2->GetID() )))
+        {
+          first  = n1->GetID();
+          second = n2->GetID();
+        }
+        else
+        {
+          myOri  = -1;
+          first  = n2->GetID();
+          second = n1->GetID();
+        }
+        myOri *= ori;
+      }
+      else
+      {
+        myOri = first = second = 0;
+      }
+    }
+    //int Node1() const { return myOri == -1 ? second : first; }
+
+    //bool IsSameOri( const std::pair<int,int>& link ) const { return link.first == Node1(); }
+  };
+}
+
 //=======================================================================
 //function : IsFaceExternal
 //purpose  : Check normal orientation of a given face
@@ -971,39 +995,179 @@ bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) const
   if ( myExternalFaces || !myVolume )
     return true;
 
-  if (myVolume->IsPoly()) {
-    XYZ aNormal, baryCenter, p0 (myPolyedre->GetFaceNode(faceIndex + 1, 1));
-    GetFaceNormal(faceIndex, aNormal.x, aNormal.y, aNormal.z);
-    GetBaryCenter(baryCenter.x, baryCenter.y, baryCenter.z);
-    XYZ insideVec (baryCenter - p0);
-    if (insideVec.Dot(aNormal) > 0)
-      return false;
+  if ( !myPolyedre ) // all classical volumes have external facet normals
     return true;
+
+  SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
+
+  if ( myPolyFacetOri[ faceIndex ])
+    return myPolyFacetOri[ faceIndex ] > 0;
+
+  int ori = 0; // -1-in, +1-out, 0-undef
+  double minProj, maxProj;
+  if ( projectNodesToNormal( faceIndex, minProj, maxProj ))
+  {
+    // all nodes are on the same side of the facet
+    ori = ( minProj < 0 ? +1 : -1 );
+    me->myPolyFacetOri[ faceIndex ] = ori;
+
+    if ( !me->myFwdLinks.empty() ) // concave polyhedron; collect oriented links
+      for ( int i = 0; i < myCurFace.myNbNodes; ++i )
+      {
+        NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
+        me->myFwdLinks.insert( make_pair( link, link.myOri ));
+      }
+    return ori > 0;
   }
 
-  // switch ( myVolumeNbNodes ) {
-  // case 4:
-  // case 5:
-  // case 10:
-  // case 13:
-  //   // only the bottom of a reversed tetrahedron can be internal
-  //   return ( myVolForward || faceIndex != 0 );
-  // case 6:
-  // case 15:
-  // case 12:
-  //   // in a forward prism, the top is internal, in a reversed one - bottom
-  //   return ( myVolForward ? faceIndex != 1 : faceIndex != 0 );
-  // case 8:
-  // case 20:
-  // case 27: {
-  //   // in a forward hexahedron, even face normal is external, odd - internal
-  //   bool odd = faceIndex % 2;
-  //   return ( myVolForward ? !odd : odd );
+  SaveFacet savedFacet( myCurFace );
+
+  // concave polyhedron
+
+  if ( me->myFwdLinks.empty() ) // get links of the least ambiguously oriented facet
+  {
+    for ( size_t i = 0; i < myPolyFacetOri.size() && !ori; ++i )
+      ori = me->myPolyFacetOri[ i ];
+
+    if ( !ori ) // none facet is oriented yet
+    {
+      // find the least ambiguously oriented facet
+      int faceMostConvex = -1;
+      std::map< double, int > convexity2face;
+      for ( size_t iF = 0; iF < myPolyFacetOri.size() && faceMostConvex < 0; ++iF )
+      {
+        if ( projectNodesToNormal( iF, minProj, maxProj ))
+        {
+          // all nodes are on the same side of the facet
+          me->myPolyFacetOri[ iF ] = ( minProj < 0 ? +1 : -1 );
+          faceMostConvex = iF;
+        }
+        else
+        {
+          ori = ( -minProj < maxProj ? -1 : +1 );
+          double convexity = std::min( -minProj, maxProj ) / std::max( -minProj, maxProj );
+          convexity2face.insert( make_pair( convexity, iF * ori ));
+        }
+      }
+      if ( faceMostConvex < 0 ) // none facet has nodes on the same side
+      {
+        // use the least ambiguous facet
+        faceMostConvex = convexity2face.begin()->second;
+        ori = ( faceMostConvex < 0 ? -1 : +1 );
+        faceMostConvex = std::abs( faceMostConvex );
+        me->myPolyFacetOri[ faceMostConvex ] = ori;
+      }
+    }
+    // collect links of the oriented facets in me->myFwdLinks
+    for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
+    {
+      ori = me->myPolyFacetOri[ iF ];
+      if ( !ori ) continue;
+      setFace( iF );
+      for ( int i = 0; i < myCurFace.myNbNodes; ++i )
+      {
+        NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1], ori );
+        me->myFwdLinks.insert( make_pair( link, link.myOri ));
+      }
+    }
+  }
+
+  // compare orientation of links of the facet with myFwdLinks
+  ori = 0;
+  setFace( faceIndex );
+  vector< NLink > links( myCurFace.myNbNodes ), links2;
+  for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
+  {
+    NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
+    std::map<Link, int>::iterator l2o = me->myFwdLinks.find( link );
+    if ( l2o != me->myFwdLinks.end() )
+      ori = link.myOri * l2o->second * -1;
+    links[ i ] = link;
+  }
+  while ( !ori ) // the facet has no common links with already oriented facets
+  {
+    // orient and collect links of other non-oriented facets
+    for ( size_t iF = 0; iF < myPolyFacetOri.size(); ++iF )
+    {
+      if ( me->myPolyFacetOri[ iF ] ) continue; // already oriented
+      setFace( iF );
+      links2.clear();
+      ori = 0;
+      for ( int i = 0; i < myCurFace.myNbNodes && !ori; ++i )
+      {
+        NLink link( myCurFace.myNodes[i], myCurFace.myNodes[i+1] );
+        std::map<Link, int>::iterator l2o = me->myFwdLinks.find( link );
+        if ( l2o != me->myFwdLinks.end() )
+          ori = link.myOri * l2o->second * -1;
+        links2.push_back( link );
+      }
+      if ( ori ) // one more facet oriented
+      {
+        me->myPolyFacetOri[ iF ] = ori;
+        for ( size_t i = 0; i < links2.size(); ++i )
+          me->myFwdLinks.insert( make_pair( links2[i], links2[i].myOri * ori ));
+        break;
+      }
+    }
+    if ( !ori )
+      return false; // error in algorithm: infinite loop
+
+    // try to orient the facet again
+    ori = 0;
+    for ( size_t i = 0; i < links.size() && !ori; ++i )
+    {
+      std::map<Link, int>::iterator l2o = me->myFwdLinks.find( links[i] );
+      if ( l2o != me->myFwdLinks.end() )
+        ori = links[i].myOri * l2o->second * -1;
+    }
+    me->myPolyFacetOri[ faceIndex ] = ori;
+  }
+
+  return ori > 0;
+}
+
+//=======================================================================
+//function : projectNodesToNormal
+//purpose  : compute min and max projections of all nodes to normal of a facet.
+//=======================================================================
+
+bool SMDS_VolumeTool::projectNodesToNormal( int     faceIndex,
+                                            double& minProj,
+                                            double& maxProj ) const
+{
+  minProj = std::numeric_limits<double>::max();
+  maxProj = std::numeric_limits<double>::min();
+
+  XYZ normal;
+  if ( !GetFaceNormal( faceIndex, normal.x, normal.y, normal.z ))
+    return false;
+  XYZ p0 ( myCurFace.myNodes[0] );
+  for ( size_t i = 0; i < myVolumeNodes.size(); ++i )
+  {
+    if ( std::find( myCurFace.myNodes.begin() + 1,
+                    myCurFace.myNodes.end(),
+                    myVolumeNodes[ i ] ) != myCurFace.myNodes.end() )
+      continue; // node of the faceIndex-th facet
+
+    double proj = normal.Dot( XYZ( myVolumeNodes[ i ]) - p0 );
+    if ( proj < minProj ) minProj = proj;
+    if ( proj > maxProj ) maxProj = proj;
+  }
+  const double tol = 1e-7;
+  minProj += tol;
+  maxProj -= tol;
+  bool diffSize = ( minProj * maxProj < 0 );
+  // if ( diffSize )
+  // {
+  //   minProj = -minProj;
   // }
-  // default:;
+  // else if ( minProj < 0 )
+  // {
+  //   minProj = -minProj;
+  //   maxProj = -maxProj;
   // }
-  // return false;
-  return true;
+
+  return !diffSize; // ? 0 : (minProj >= 0);
 }
 
 //=======================================================================
@@ -1016,16 +1180,16 @@ bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, doub
   if ( !setFace( faceIndex ))
     return false;
 
-  const int iQuad = ( myFaceNbNodes > 6 && !myPolyedre ) ? 2 : 1;
-  XYZ p1 ( myFaceNodes[0*iQuad] );
-  XYZ p2 ( myFaceNodes[1*iQuad] );
-  XYZ p3 ( myFaceNodes[2*iQuad] );
+  const int iQuad = ( !myPolyedre && myCurFace.myNbNodes > 6 ) ? 2 : 1;
+  XYZ p1 ( myCurFace.myNodes[0*iQuad] );
+  XYZ p2 ( myCurFace.myNodes[1*iQuad] );
+  XYZ p3 ( myCurFace.myNodes[2*iQuad] );
   XYZ aVec12( p2 - p1 );
   XYZ aVec13( p3 - p1 );
   XYZ cross = aVec12.Crossed( aVec13 );
 
-  if ( myFaceNbNodes >3*iQuad ) {
-    XYZ p4 ( myFaceNodes[3*iQuad] );
+  if ( myCurFace.myNbNodes >3*iQuad ) {
+    XYZ p4 ( myCurFace.myNodes[3*iQuad] );
     XYZ aVec14( p4 - p1 );
     XYZ cross2 = aVec13.Crossed( aVec14 );
     cross = cross + cross2;
@@ -1054,11 +1218,11 @@ bool SMDS_VolumeTool::GetFaceBaryCenter (int faceIndex, double & X, double & Y,
     return false;
 
   X = Y = Z = 0.0;
-  for ( int i = 0; i < myFaceNbNodes; ++i )
+  for ( int i = 0; i < myCurFace.myNbNodes; ++i )
   {
-    X += myFaceNodes[i]->X() / myFaceNbNodes;
-    Y += myFaceNodes[i]->Y() / myFaceNbNodes;
-    Z += myFaceNodes[i]->Z() / myFaceNbNodes;
+    X += myCurFace.myNodes[i]->X() / myCurFace.myNbNodes;
+    Y += myCurFace.myNodes[i]->Y() / myCurFace.myNbNodes;
+    Z += myCurFace.myNodes[i]->Z() / myCurFace.myNbNodes;
   }
   return true;
 }
@@ -1070,27 +1234,36 @@ bool SMDS_VolumeTool::GetFaceBaryCenter (int faceIndex, double & X, double & Y,
 
 double SMDS_VolumeTool::GetFaceArea( int faceIndex ) const
 {
-  if (myVolume->IsPoly()) {
-    MESSAGE("Warning: attempt to obtain area of a face of polyhedral volume");
-    return 0;
-  }
-
+  double area = 0;
   if ( !setFace( faceIndex ))
-    return 0;
+    return area;
 
-  XYZ p1 ( myFaceNodes[0] );
-  XYZ p2 ( myFaceNodes[1] );
-  XYZ p3 ( myFaceNodes[2] );
+  XYZ p1 ( myCurFace.myNodes[0] );
+  XYZ p2 ( myCurFace.myNodes[1] );
+  XYZ p3 ( myCurFace.myNodes[2] );
   XYZ aVec12( p2 - p1 );
   XYZ aVec13( p3 - p1 );
-  double area = aVec12.Crossed( aVec13 ).Magnitude() * 0.5;
+  area += aVec12.Crossed( aVec13 ).Magnitude();
 
-  if ( myFaceNbNodes == 4 ) {
-    XYZ p4 ( myFaceNodes[3] );
-    XYZ aVec14( p4 - p1 );
-    area += aVec14.Crossed( aVec13 ).Magnitude() * 0.5;
+  if (myVolume->IsPoly())
+  {
+    for ( int i = 3; i < myCurFace.myNbNodes; ++i )
+    {
+      XYZ pI ( myCurFace.myNodes[i] );
+      XYZ aVecI( pI - p1 );
+      area += aVec13.Crossed( aVecI ).Magnitude();
+      aVec13 = aVecI;
+    }
   }
-  return area;
+  else
+  {
+    if ( myCurFace.myNbNodes == 4 ) {
+      XYZ p4 ( myCurFace.myNodes[3] );
+      XYZ aVec14( p4 - p1 );
+      area += aVec14.Crossed( aVec13 ).Magnitude();
+    }
+  }
+  return area / 2;
 }
 
 //================================================================================
@@ -1101,7 +1274,7 @@ double SMDS_VolumeTool::GetFaceArea( int faceIndex ) const
 
 int SMDS_VolumeTool::GetCenterNodeIndex( int faceIndex ) const
 {
-  if ( myAllFacesNbNodes && myVolumeNbNodes == 27 ) // classic element with 27 nodes
+  if ( myAllFacesNbNodes && myVolumeNodes.size() == 27 ) // classic element with 27 nodes
   {
     switch ( faceIndex ) {
     case 0: return 20;
@@ -1129,7 +1302,7 @@ int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
   const int nbHoriFaces = 2;
 
   if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
-    switch ( myVolumeNbNodes ) {
+    switch ( myVolumeNodes.size() ) {
     case 6:
     case 15:
       if ( faceIndex == 0 || faceIndex == 1 )
@@ -1182,31 +1355,42 @@ bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
       MESSAGE("Warning: bad volumic element");
       return false;
     }
-    bool isLinked = false;
-    int iface;
-    for (iface = 1; iface <= myNbFaces && !isLinked; iface++) {
-      int inode, nbFaceNodes = myPolyedre->NbFaceNodes(iface);
-
-      for (inode = 1; inode <= nbFaceNodes && !isLinked; inode++) {
-        const SMDS_MeshNode* curNode = myPolyedre->GetFaceNode(iface, inode);
-
-        if (curNode == theNode1 || curNode == theNode2) {
-          int inextnode = (inode == nbFaceNodes) ? 1 : inode + 1;
-          const SMDS_MeshNode* nextNode = myPolyedre->GetFaceNode(iface, inextnode);
-
-          if ((curNode == theNode1 && nextNode == theNode2) ||
-              (curNode == theNode2 && nextNode == theNode1)) {
-            isLinked = true;
-          }
-        }
+    if ( !myAllFacesNbNodes ) {
+      SMDS_VolumeTool*  me = const_cast< SMDS_VolumeTool* >( this );
+      me->myPolyQuantities = myPolyedre->GetQuantities();
+      myAllFacesNbNodes    = &myPolyQuantities[0];
+    }
+    int from, to = 0, d1 = 1, d2 = 2;
+    if ( myPolyedre->IsQuadratic() ) {
+      if ( theIgnoreMediumNodes ) {
+        d1 = 2; d2 = 0;
+      }
+    } else {
+      d2 = 0;
+    }
+    vector<const SMDS_MeshNode*>::const_iterator i;
+    for (int iface = 0; iface < myNbFaces; iface++)
+    {
+      from = to;
+      to  += myPolyQuantities[iface];
+      i    = std::find( myVolumeNodes.begin() + from, myVolumeNodes.begin() + to, theNode1 );
+      if ( i != myVolumeNodes.end() )
+      {
+        if ((  theNode2 == *( i-d1 ) ||
+               theNode2 == *( i+d1 )))
+          return true;
+        if (( d2 ) &&
+            (( theNode2 == *( i-d2 ) ||
+               theNode2 == *( i+d2 ))))
+          return true;
       }
     }
-    return isLinked;
+    return false;
   }
 
   // find nodes indices
   int i1 = -1, i2 = -1, nbFound = 0;
-  for ( int i = 0; i < myVolumeNbNodes && nbFound < 2; i++ )
+  for ( int i = 0; i < myVolumeNodes.size() && nbFound < 2; i++ )
   {
     if ( myVolumeNodes[ i ] == theNode1 )
       i1 = i, ++nbFound;
@@ -1234,7 +1418,7 @@ bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
   int minInd = min( theNode1Index, theNode2Index );
   int maxInd = max( theNode1Index, theNode2Index );
 
-  if ( minInd < 0 || maxInd > myVolumeNbNodes - 1 || maxInd == minInd )
+  if ( minInd < 0 || maxInd > myVolumeNodes.size() - 1 || maxInd == minInd )
     return false;
 
   VolumeType type = GetVolumeType();
@@ -1351,7 +1535,7 @@ bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
 {
   if ( myVolume ) {
-    for ( int i = 0; i < myVolumeNbNodes; i++ ) {
+    for ( int i = 0; i < myVolumeNodes.size(); i++ ) {
       if ( myVolumeNodes[ i ] == theNode )
         return i;
     }
@@ -1370,24 +1554,32 @@ int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
 int SMDS_VolumeTool::GetAllExistingFaces(vector<const SMDS_MeshElement*> & faces) const
 {
   faces.clear();
-  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;
+  SaveFacet savedFacet( myCurFace );
+  if ( IsPoly() )
+    for ( int iF = 0; iF < NbFaces(); ++iF ) {
+      if ( setFace( iF ))
+        if ( const SMDS_MeshElement* face = SMDS_Mesh::FindFace( myCurFace.myNodes ))
+          faces.push_back( face );
+    }
+  else
+    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 );
     }
-    if ( face )
-      faces.push_back( face );
-  }
   return faces.size();
 }
 
@@ -1395,17 +1587,17 @@ int SMDS_VolumeTool::GetAllExistingFaces(vector<const SMDS_MeshElement*> & faces
 //================================================================================
 /*!
  * \brief Fill vector with boundary edges existing in the mesh
-  * \param edges - vector of found edges
-  * \retval int - nb of found faces
+ * \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-1; ++i ) {
-    for ( int j = i + 1; j < myVolumeNbNodes; ++j ) {
+  edges.reserve( myVolumeNodes.size() * 2 );
+  for ( int i = 0; i < myVolumeNodes.size()-1; ++i ) {
+    for ( int j = i + 1; j < myVolumeNodes.size(); ++j ) {
       if ( IsLinked( i, j )) {
         const SMDS_MeshElement* edge =
           SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] );
@@ -1428,30 +1620,20 @@ 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;
-  
+  SaveFacet savedFacet( myCurFace );
+
   // it seems that compute distance twice is faster than organization of a sole computing
-  myCurFace = -1;
+  myCurFace.myIndex = -1;
   for ( int iF = 0; iF < myNbFaces; ++iF )
   {
     setFace( iF );
-    for ( int iN = 0; iN < myFaceNbNodes; iN += iQ )
+    for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
     {
-      XYZ n1( myFaceNodes[ iN ]);
-      XYZ n2( myFaceNodes[(iN + iQ) % myFaceNbNodes]);
+      XYZ n1( myCurFace.myNodes[ iN ]);
+      XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
       minSize = std::min( minSize, (n1 - n2).SquareMagnitude());
     }
   }
-  // restore current face data
-  myCurFace = curFace;
-  myFaceNbNodes = nbN;
-  myFaceNodeIndices = ind;
-  delete [] myFaceNodes; myFaceNodes = nodes;
 
   return minSize;
 }
@@ -1467,30 +1649,20 @@ double SMDS_VolumeTool::MaxLinearSize2() const
   double maxSize = -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;
+  SaveFacet savedFacet( myCurFace );
   
   // it seems that compute distance twice is faster than organization of a sole computing
-  myCurFace = -1;
+  myCurFace.myIndex = -1;
   for ( int iF = 0; iF < myNbFaces; ++iF )
   {
     setFace( iF );
-    for ( int iN = 0; iN < myFaceNbNodes; iN += iQ )
+    for ( int iN = 0; iN < myCurFace.myNbNodes; iN += iQ )
     {
-      XYZ n1( myFaceNodes[ iN ]);
-      XYZ n2( myFaceNodes[(iN + iQ) % myFaceNbNodes]);
+      XYZ n1( myCurFace.myNodes[ iN ]);
+      XYZ n2( myCurFace.myNodes[(iN + iQ) % myCurFace.myNbNodes]);
       maxSize = std::max( maxSize, (n1 - n2).SquareMagnitude());
     }
   }
-  // restore current face data
-  myCurFace         = curFace;
-  myFaceNbNodes     = nbN;
-  myFaceNodeIndices = ind;
-  delete [] myFaceNodes; myFaceNodes = nodes;
 
   return maxSize;
 }
@@ -1506,13 +1678,13 @@ bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherV
 {
   const bool isFree = true;
 
-  if (!setFace( faceIndex ))
+  if ( !setFace( faceIndex ))
     return !isFree;
 
   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
 
   const int  di = myVolume->IsQuadratic() ? 2 : 1;
-  const int nbN = ( myFaceNbNodes/di <= 4 && !IsPoly()) ? 3 : myFaceNbNodes/di; // nb nodes to check
+  const int nbN = ( myCurFace.myNbNodes/di <= 4 && !IsPoly()) ? 3 : myCurFace.myNbNodes/di; // nb nodes to check
 
   SMDS_ElemIteratorPtr eIt = nodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
   while ( eIt->more() )
@@ -1528,7 +1700,7 @@ bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherV
     {
       // if ( vol->IsPoly() || vol->NbFaces() > 6 ) // vol is polyhed or hex prism 
       // {
-      //   int nb = myFaceNbNodes;
+      //   int nb = myCurFace.myNbNodes;
       //   if ( myVolume->GetEntityType() != vol->GetEntityType() )
       //     nb -= ( GetCenterNodeIndex(0) > 0 );
       //   set<const SMDS_MeshNode*> faceNodes( nodes, nodes + nb );
@@ -1557,7 +1729,7 @@ bool SMDS_VolumeTool::IsFreeFaceAdv( int faceIndex, const SMDS_MeshElement** oth
     return !isFree;
 
   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
-  const int nbFaceNodes = myFaceNbNodes;
+  const int nbFaceNodes = myCurFace.myNbNodes;
 
   // evaluate nb of face nodes shared by other volumes
   int maxNbShared = -1;
@@ -1736,19 +1908,14 @@ bool SMDS_VolumeTool::setFace( int faceIndex ) const
   if ( !myVolume )
     return false;
 
-  if ( myCurFace == faceIndex )
+  if ( myCurFace.myIndex == faceIndex )
     return true;
 
-  myCurFace = -1;
+  myCurFace.myIndex = -1;
 
   if ( faceIndex < 0 || faceIndex >= NbFaces() )
     return false;
 
-  if (myFaceNodes != NULL) {
-    delete [] myFaceNodes;
-    myFaceNodes = NULL;
-  }
-
   if (myVolume->IsPoly())
   {
     if (!myPolyedre) {
@@ -1757,21 +1924,31 @@ bool SMDS_VolumeTool::setFace( int faceIndex ) const
     }
 
     // set face nodes
-    int iNode;
-    myFaceNbNodes = myPolyedre->NbFaceNodes(faceIndex + 1);
-    myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
-    for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
-      myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, iNode + 1);
-    myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ]; // last = first
+    SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* >( this );
+    if ( !myAllFacesNbNodes ) {
+      me->myPolyQuantities = myPolyedre->GetQuantities();
+      myAllFacesNbNodes    = &myPolyQuantities[0];
+    }
+    myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
+    myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
+    me->myPolyIndices.resize( myCurFace.myNbNodes + 1 );
+    myCurFace.myNodeIndices = & me->myPolyIndices[0];
+    int shift = std::accumulate( myAllFacesNbNodes, myAllFacesNbNodes+faceIndex, 0 );
+    for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
+    {
+      myCurFace.myNodes      [ iNode ] = myVolumeNodes[ shift + iNode ];
+      myCurFace.myNodeIndices[ iNode ] = shift + iNode;
+    }
+    myCurFace.myNodes      [ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ]; // last = first
+    myCurFace.myNodeIndices[ myCurFace.myNbNodes ] = myCurFace.myNodeIndices[ 0 ];
 
     // check orientation
     if (myExternalFaces)
     {
-      myCurFace = faceIndex; // avoid infinite recursion in IsFaceExternal()
+      myCurFace.myIndex = 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] );
+        std::reverse( myCurFace.myNodes.begin(), myCurFace.myNodes.end() );
       myExternalFaces = true;
     }
   }
@@ -1780,7 +1957,7 @@ bool SMDS_VolumeTool::setFace( int faceIndex ) const
     if ( !myAllFacesNodeIndices_F )
     {
       // choose data for an element type
-      switch ( myVolumeNbNodes ) {
+      switch ( myVolumeNodes.size() ) {
       case 4:
         myAllFacesNodeIndices_F  = &Tetra_F [0][0];
         //myAllFacesNodeIndices_FE = &Tetra_F [0][0];
@@ -1837,7 +2014,7 @@ bool SMDS_VolumeTool::setFace( int faceIndex ) const
         myAllFacesNodeIndices_RE = &QuadHexa_RE[0][0];
         myAllFacesNbNodes        = QuadHexa_nbN;
         myMaxFaceNbNodes         = sizeof(QuadHexa_F[0])/sizeof(QuadHexa_F[0][0]);
-        if ( !myIgnoreCentralNodes && myVolumeNbNodes == 27 )
+        if ( !myIgnoreCentralNodes && myVolumeNodes.size() == 27 )
         {
           myAllFacesNodeIndices_F  = &TriQuadHexa_F [0][0];
           //myAllFacesNodeIndices_FE = &TriQuadHexa_FE[0][0];
@@ -1857,21 +2034,21 @@ bool SMDS_VolumeTool::setFace( int faceIndex ) const
         return false;
       }
     }
-    myFaceNbNodes = myAllFacesNbNodes[ faceIndex ];
+    myCurFace.myNbNodes = myAllFacesNbNodes[ faceIndex ];
     // if ( myExternalFaces )
-    //   myFaceNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
+    //   myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
     // else
-    //   myFaceNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes );
-    myFaceNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
+    //   myCurFace.myNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes );
+    myCurFace.myNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
 
     // set face nodes
-    myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
-    for ( int iNode = 0; iNode < myFaceNbNodes; iNode++ )
-      myFaceNodes[ iNode ] = myVolumeNodes[ myFaceNodeIndices[ iNode ]];
-    myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ];
+    myCurFace.myNodes.resize( myCurFace.myNbNodes + 1 );
+    for ( int iNode = 0; iNode < myCurFace.myNbNodes; iNode++ )
+      myCurFace.myNodes[ iNode ] = myVolumeNodes[ myCurFace.myNodeIndices[ iNode ]];
+    myCurFace.myNodes[ myCurFace.myNbNodes ] = myCurFace.myNodes[ 0 ];
   }
 
-  myCurFace = faceIndex;
+  myCurFace.myIndex = faceIndex;
 
   return true;
 }
diff --git a/src/SMDS/SMDS_VolumeTool.hxx b/src/SMDS/SMDS_VolumeTool.hxx
index 4b0d63e2f..7e88f7c1a 100644
--- a/src/SMDS/SMDS_VolumeTool.hxx
+++ b/src/SMDS/SMDS_VolumeTool.hxx
@@ -38,6 +38,7 @@ class SMDS_MeshVolume;
 
 #include <vector>
 #include <set>
+#include <map>
 
 // =========================================================================
 //
@@ -90,10 +91,10 @@ class SMDS_EXPORT SMDS_VolumeTool
   // top and bottom faces are reversed.
   // Result of IsForward() and methods returning nodes change
 
-  const SMDS_MeshNode** GetNodes() { return myVolumeNodes; }
+  const SMDS_MeshNode** GetNodes() { return &myVolumeNodes[0]; }
   // Return array of volume nodes
 
-  int NbNodes() { return myVolumeNbNodes; }
+  int NbNodes() { return myVolumeNodes.size(); }
   // Return array of volume nodes
 
   double GetSize() const;
@@ -242,31 +243,42 @@ class SMDS_EXPORT SMDS_VolumeTool
   static int GetOppFaceIndexOfHex( int faceIndex );
   // Return index of the opposite face of the hexahedron
 
-private:
+ private:
 
   bool setFace( int faceIndex ) const;
 
+  bool projectNodesToNormal( int faceIndex, double& minProj, double& maxProj ) const;
+
   const SMDS_MeshElement* myVolume;
   const SMDS_VtkVolume*   myPolyedre;
   bool                    myIgnoreCentralNodes;
 
   bool                    myVolForward;
   int                     myNbFaces;
-  int                     myVolumeNbNodes;
-  const SMDS_MeshNode**   myVolumeNodes;
-  std::vector< int >      myPolyIndices;
+  std::vector<const SMDS_MeshNode*> myVolumeNodes;
+  std::vector< int >      myPolyIndices; // of a myCurFace
+  std::vector< int >      myPolyQuantities;
+  std::vector< int >      myPolyFacetOri; // -1-in, +1-out, 0-undef
 
-  mutable bool                    myExternalFaces;
+  typedef std::pair<int,int> Link;
+  std::map<Link, int>     myFwdLinks; // used in IsFaceExternal() to find out myPolyFacetOri
 
-  mutable const int*              myAllFacesNodeIndices_F;
-  mutable const int*              myAllFacesNodeIndices_RE;
-  mutable const int*              myAllFacesNbNodes;
-  mutable int                     myMaxFaceNbNodes;
+  mutable bool            myExternalFaces;
 
-  mutable int                     myCurFace;
-  mutable int                     myFaceNbNodes;
-  mutable int*                    myFaceNodeIndices;
-  mutable const SMDS_MeshNode**   myFaceNodes;
+  mutable const int*      myAllFacesNodeIndices_F;
+  mutable const int*      myAllFacesNodeIndices_RE;
+  mutable const int*      myAllFacesNbNodes;
+  mutable int             myMaxFaceNbNodes;
+
+  struct SaveFacet;
+  struct Facet
+  {
+    int                   myIndex;
+    int                   myNbNodes;
+    int*                  myNodeIndices;
+    std::vector<const SMDS_MeshNode*> myNodes;
+  };
+  mutable Facet           myCurFace;
 
 };
 #endif
diff --git a/src/SMESH/SMESH_MeshEditor.cxx b/src/SMESH/SMESH_MeshEditor.cxx
index 9b8ead42d..44bd3e24c 100644
--- a/src/SMESH/SMESH_MeshEditor.cxx
+++ b/src/SMESH/SMESH_MeshEditor.cxx
@@ -12577,6 +12577,11 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
   TConnectivity tgtNodes;
   ElemFeatures elemKind( missType ), elemToCopy;
 
+  vector<const SMDS_MeshElement*> presentBndElems;
+  vector<TConnectivity>           missingBndElems;
+  vector<int>                     freeFacets;
+  TConnectivity nodes, elemNodes;
+
   SMDS_ElemIteratorPtr eIt;
   if (elements.empty()) eIt = aMesh->elementsIterator(elemType);
   else                  eIt = elemSetIterator( elements );
@@ -12590,18 +12595,24 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
     // ------------------------------------------------------------------------------------
     // 1. For an elem, get present bnd elements and connectivities of missing bnd elements
     // ------------------------------------------------------------------------------------
-    vector<const SMDS_MeshElement*> presentBndElems;
-    vector<TConnectivity>           missingBndElems;
-    TConnectivity nodes, elemNodes;
+    presentBndElems.clear();
+    missingBndElems.clear();
+    freeFacets.clear(); nodes.clear(); elemNodes.clear();
     if ( vTool.Set(elem, /*ignoreCentralNodes=*/true) ) // elem is a volume --------------
     {
-      vTool.SetExternalNormal();
       const SMDS_MeshElement* otherVol = 0;
       for ( int iface = 0, n = vTool.NbFaces(); iface < n; iface++ )
       {
         if ( !vTool.IsFreeFace(iface, &otherVol) &&
              ( !aroundElements || elements.count( otherVol )))
           continue;
+        freeFacets.push_back( iface );
+      }
+      if ( missType == SMDSAbs_Face )
+        vTool.SetExternalNormal();
+      for ( size_t i = 0; i < freeFacets.size(); ++i )
+      {
+        int                iface = freeFacets[i];
         const SMDS_MeshNode** nn = vTool.GetFaceNodes(iface);
         const size_t nbFaceNodes = vTool.NbFaceNodes (iface);
         if ( missType == SMDSAbs_Edge ) // boundary edges
diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx
index aec095ec3..ab24f4c20 100644
--- a/src/SMESH/SMESH_MesherHelper.cxx
+++ b/src/SMESH/SMESH_MesherHelper.cxx
@@ -807,6 +807,24 @@ GeomAPI_ProjectPointOnSurf& SMESH_MesherHelper::GetProjector(const TopoDS_Face&
   return *( i_proj->second );
 }
 
+//=======================================================================
+//function : GetSurface
+//purpose  : Return a cached ShapeAnalysis_Surface of a FACE
+//=======================================================================
+
+Handle(ShapeAnalysis_Surface) SMESH_MesherHelper::GetSurface(const TopoDS_Face& F )
+{
+  Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
+  int faceID = GetMeshDS()->ShapeToIndex( F );
+  TID2Surface::iterator i_surf = myFace2Surface.find( faceID );
+  if ( i_surf == myFace2Surface.end() && faceID )
+  {
+    Handle(ShapeAnalysis_Surface) surf( new ShapeAnalysis_Surface( surface ));
+    i_surf = myFace2Surface.insert( make_pair( faceID, surf )).first;
+  }
+  return i_surf->second;
+}
+
 namespace
 {
   gp_XY AverageUV(const gp_XY& uv1, const gp_XY& uv2) { return ( uv1 + uv2 ) / 2.; }
@@ -1581,6 +1599,26 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
   {
     F = TopoDS::Face(meshDS->IndexToShape( faceID = pos.first ));
     uv[0] = GetNodeUV(F,n1,n2, force3d ? 0 : &uvOK[0]);
+    if ( HasDegeneratedEdges() && !force3d ) // IPAL52850 (degen VERTEX not at singularity)
+    {
+      // project middle point to a surface
+      SMESH_TNodeXYZ p1( n1 ), p2( n2 );
+      gp_Pnt pMid = 0.5 * ( p1 + p2 );
+      Handle(ShapeAnalysis_Surface) projector = GetSurface( F );
+      gp_Pnt2d uvMid;
+      if ( uvOK[0] )
+        uvMid = projector->NextValueOfUV( uv[0], pMid, BRep_Tool::Tolerance( F ));
+      else
+        uvMid = projector->ValueOfUV( pMid, getFaceMaxTol( F ));
+      if ( projector->Gap() * projector->Gap() < ( p1 - p2 ).SquareModulus() / 4 )
+      {
+        gp_Pnt pProj = projector->Value( uvMid );
+        n12  = meshDS->AddNode( pProj.X(), pProj.Y(), pProj.Z() );
+        meshDS->SetNodeOnFace( n12, faceID, uvMid.X(), uvMid.Y() );
+        myTLinkNodeMap.insert( make_pair ( link, n12 ));
+        return n12;
+      }
+    }
     uv[1] = GetNodeUV(F,n2,n1, force3d ? 0 : &uvOK[1]);
   }
   else if ( pos.second == TopAbs_EDGE )
diff --git a/src/SMESH/SMESH_MesherHelper.hxx b/src/SMESH/SMESH_MesherHelper.hxx
index d85f8329e..9e8d82958 100644
--- a/src/SMESH/SMESH_MesherHelper.hxx
+++ b/src/SMESH/SMESH_MesherHelper.hxx
@@ -34,6 +34,7 @@
 #include <SMDS_QuadraticEdge.hxx>
 
 #include <Geom_Surface.hxx>
+#include <ShapeAnalysis_Surface.hxx>
 #include <TopoDS_Face.hxx>
 #include <TopoDS_Shape.hxx>
 #include <gp_Pnt2d.hxx>
@@ -528,6 +529,10 @@ public:
   GeomAPI_ProjectPointOnSurf& GetProjector(const TopoDS_Face& F,
                                            TopLoc_Location&   loc,
                                            double             tol=0 ) const; 
+  /*!
+   * \brief Return a cached ShapeAnalysis_Surface of a FACE
+   */
+  Handle(ShapeAnalysis_Surface) GetSurface(const TopoDS_Face& F ); 
 
   /*!
    * \brief Check if shape is a degenerated edge or it's vertex
@@ -731,8 +736,10 @@ public:
 
   std::map< int, double > myFaceMaxTol;
 
+  typedef std::map< int, Handle(ShapeAnalysis_Surface)> TID2Surface;
   typedef std::map< int, GeomAPI_ProjectPointOnSurf* >  TID2ProjectorOnSurf;
   typedef std::map< int, GeomAPI_ProjectPointOnCurve* > TID2ProjectorOnCurve;
+  TID2Surface          myFace2Surface;
   TID2ProjectorOnSurf  myFace2Projector;
   TID2ProjectorOnCurve myEdge2Projector;
 
diff --git a/src/SMESH_SWIG/smeshBuilder.py b/src/SMESH_SWIG/smeshBuilder.py
index 3810ea576..6cf59af32 100644
--- a/src/SMESH_SWIG/smeshBuilder.py
+++ b/src/SMESH_SWIG/smeshBuilder.py
@@ -3619,8 +3619,8 @@ class Mesh:
     ## Creates 2D mesh as skin on boundary faces of a 3D mesh
     #  @return TRUE if operation has been completed successfully, FALSE otherwise
     #  @ingroup l2_modif_edit
-    def  Make2DMeshFrom3D(self):
-        return self.editor. Make2DMeshFrom3D()
+    def Make2DMeshFrom3D(self):
+        return self.editor.Make2DMeshFrom3D()
 
     ## Creates missing boundary elements
     #  @param elements - elements whose boundary is to be checked:
@@ -4910,6 +4910,8 @@ class Mesh:
         fun = self._getFunctor( funType )
         if fun:
             if meshPart:
+                if hasattr( meshPart, "SetMesh" ):
+                    meshPart.SetMesh( self.mesh ) # set mesh to filter
                 hist = fun.GetLocalHistogram( 1, False, meshPart )
             else:
                 hist = fun.GetHistogram( 1, False )
diff --git a/src/StdMeshers/StdMeshers_FaceSide.cxx b/src/StdMeshers/StdMeshers_FaceSide.cxx
index 5882b13af..2625b9fb1 100644
--- a/src/StdMeshers/StdMeshers_FaceSide.cxx
+++ b/src/StdMeshers/StdMeshers_FaceSide.cxx
@@ -410,7 +410,7 @@ const vector<UVPtStruct>& StdMeshers_FaceSide::GetUVPtStruct(bool   isXConst,
       if ( u2node.empty() ) return myPoints;
 
       const SMDS_MeshNode* node;
-      if ( IsClosed() )
+      if ( IsClosed() && !proxySubMesh[0] )
         node = u2node.begin()->second;
       else
       {