Salome HOME
Regression revealed after restoring mesh data destruction at study closing
[modules/smesh.git] / src / SMESH / SMESH_MeshEditor.cxx
index 56371d7112930ca2ba49711fe92bed76519cbf84..ceaeb1acb681d86934f3c211e4b07704a6a7622b 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2016  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
 
 #include "SMESH_MeshEditor.hxx"
 
-#include "SMDS_FaceOfNodes.hxx"
-#include "SMDS_VolumeTool.hxx"
+#include "SMDS_Downward.hxx"
 #include "SMDS_EdgePosition.hxx"
+#include "SMDS_FaceOfNodes.hxx"
 #include "SMDS_FacePosition.hxx"
-#include "SMDS_SpacePosition.hxx"
-#include "SMDS_MeshGroup.hxx"
 #include "SMDS_LinearEdge.hxx"
-#include "SMDS_Downward.hxx"
+#include "SMDS_MeshGroup.hxx"
 #include "SMDS_SetIterator.hxx"
-
+#include "SMDS_SpacePosition.hxx"
+#include "SMDS_VolumeTool.hxx"
 #include "SMESHDS_Group.hxx"
 #include "SMESHDS_Mesh.hxx"
-
 #include "SMESH_Algo.hxx"
 #include "SMESH_ControlsDef.hxx"
 #include "SMESH_Group.hxx"
+#include "SMESH_Mesh.hxx"
 #include "SMESH_MeshAlgos.hxx"
 #include "SMESH_MesherHelper.hxx"
 #include "SMESH_OctreeNode.hxx"
 #include "SMESH_subMesh.hxx"
 
-#include <Basics_OCCTVersion.hxx>
-
 #include "utilities.h"
 #include "chrono.hxx"
 
@@ -73,6 +70,7 @@
 #include <TopTools_ListOfShape.hxx>
 #include <TopTools_SequenceOfShape.hxx>
 #include <TopoDS.hxx>
+#include <TopoDS_Edge.hxx>
 #include <TopoDS_Face.hxx>
 #include <TopoDS_Solid.hxx>
 #include <gp.hxx>
 #include <sstream>
 
 #include <boost/tuple/tuple.hpp>
+#include <boost/container/flat_set.hpp>
 
 #include <Standard_Failure.hxx>
 #include <Standard_ErrorHandler.hxx>
+#include <OSD_Parallel.hxx>
+
+#include "SMESH_TryCatch.hxx" // include after OCCT headers!
 
 #define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
 
 using namespace std;
 using namespace SMESH::Controls;
 
-namespace
-{
-  template < class ELEM_SET >
-  SMDS_ElemIteratorPtr elemSetIterator( const ELEM_SET& elements )
-  {
-    typedef SMDS_SetIterator
-      < SMDS_pElement, typename ELEM_SET::const_iterator> TSetIterator;
-    return SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() ));
-  }
-}
-
 //=======================================================================
 //function : SMESH_MeshEditor
 //purpose  :
@@ -125,6 +116,18 @@ SMESH_MeshEditor::SMESH_MeshEditor( SMESH_Mesh* theMesh )
 {
 }
 
+//================================================================================
+/*!
+ * \brief Return mesh DS
+ */
+//================================================================================
+
+SMESHDS_Mesh * SMESH_MeshEditor::GetMeshDS()
+{
+  return myMesh->GetMeshDS();
+}
+
+
 //================================================================================
 /*!
  * \brief Clears myLastCreatedNodes and myLastCreatedElems
@@ -133,8 +136,8 @@ SMESH_MeshEditor::SMESH_MeshEditor( SMESH_Mesh* theMesh )
 
 void SMESH_MeshEditor::ClearLastCreated()
 {
-  myLastCreatedNodes.Clear();
-  myLastCreatedElems.Clear();
+  SMESHUtils::FreeVector( myLastCreatedElems );
+  SMESHUtils::FreeVector( myLastCreatedNodes );
 }
 
 //================================================================================
@@ -159,7 +162,7 @@ SMESH_MeshEditor::ElemFeatures::Init( const SMDS_MeshElement* elem, bool basicOn
         myIsQuad = elem->IsQuadratic();
         if ( myType == SMDSAbs_Volume && !basicOnly )
         {
-          vector<int > quant = static_cast<const SMDS_VtkVolume* >( elem )->GetQuantities();
+          vector<int> quant = static_cast<const SMDS_MeshVolume* >( elem )->GetQuantities();
           myPolyhedQuantities.swap( quant );
         }
       }
@@ -366,7 +369,7 @@ SMESH_MeshEditor::AddElement(const vector<const SMDS_MeshNode*> & node,
 
   default:;
   }
-  if ( e ) myLastCreatedElems.Append( e );
+  if ( e ) myLastCreatedElems.push_back( e );
   return e;
 }
 
@@ -400,8 +403,7 @@ SMDS_MeshElement* SMESH_MeshEditor::AddElement(const vector<int> & nodeIDs,
 int SMESH_MeshEditor::Remove (const list< int >& theIDs,
                               const bool         isNodes )
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
   SMESHDS_Mesh* aMesh = GetMeshDS();
   set< SMESH_subMesh *> smmap;
@@ -462,31 +464,26 @@ int SMESH_MeshEditor::Remove (const list< int >& theIDs,
 
 //================================================================================
 /*!
- * \brief Create 0D elements on all nodes of the given object except those
- *        nodes on which a 0D element already exists.
+ * \brief Create 0D elements on all nodes of the given object.
  *  \param elements - Elements on whose nodes to create 0D elements; if empty, 
  *                    the all mesh is treated
  *  \param all0DElems - returns all 0D elements found or created on nodes of \a elements
+ *  \param duplicateElements - to add one more 0D element to a node or not
  */
 //================================================================================
 
 void SMESH_MeshEditor::Create0DElementsOnAllNodes( const TIDSortedElemSet& elements,
-                                                   TIDSortedElemSet&       all0DElems )
+                                                   TIDSortedElemSet&       all0DElems,
+                                                   const bool              duplicateElements )
 {
   SMDS_ElemIteratorPtr elemIt;
-  vector< const SMDS_MeshElement* > allNodes;
   if ( elements.empty() )
   {
-    allNodes.reserve( GetMeshDS()->NbNodes() );
     elemIt = GetMeshDS()->elementsIterator( SMDSAbs_Node );
-    while ( elemIt->more() )
-      allNodes.push_back( elemIt->next() );
-
-    elemIt = elemSetIterator( allNodes );
   }
   else
   {
-    elemIt = elemSetIterator( elements );
+    elemIt = SMESHUtils::elemSetIterator( elements );
   }
 
   while ( elemIt->more() )
@@ -497,12 +494,13 @@ void SMESH_MeshEditor::Create0DElementsOnAllNodes( const TIDSortedElemSet& eleme
     {
       const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
       SMDS_ElemIteratorPtr it0D = n->GetInverseElementIterator( SMDSAbs_0DElement );
-      if ( it0D->more() )
-        all0DElems.insert( it0D->next() );
-      else {
-        myLastCreatedElems.Append( GetMeshDS()->Add0DElement( n ));
-        all0DElems.insert( myLastCreatedElems.Last() );
+      if ( duplicateElements || !it0D->more() )
+      {
+        myLastCreatedElems.push_back( GetMeshDS()->Add0DElement( n ));
+        all0DElems.insert( myLastCreatedElems.back() );
       }
+      while ( it0D->more() )
+        all0DElems.insert( it0D->next() );
     }
   }
 }
@@ -515,8 +513,7 @@ void SMESH_MeshEditor::Create0DElementsOnAllNodes( const TIDSortedElemSet& eleme
 
 int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
   SMESHDS_Mesh * aMesh = GetMeshDS();
   if ( aMesh->ShapeToMesh().IsNull() )
@@ -683,20 +680,18 @@ static bool getNodesFromTwoTria(const SMDS_MeshElement * theTria1,
 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1,
                                     const SMDS_MeshElement * theTria2 )
 {
-  MESSAGE("InverseDiag");
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
-  if (!theTria1 || !theTria2)
+  if ( !theTria1 || !theTria2 ||
+       !dynamic_cast<const SMDS_MeshCell*>( theTria1 ) ||
+       !dynamic_cast<const SMDS_MeshCell*>( theTria2 ) ||
+       theTria1->GetType() != SMDSAbs_Face ||
+       theTria2->GetType() != SMDSAbs_Face )
     return false;
 
-  const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( theTria1 );
-  if (!F1) return false;
-  const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( theTria2 );
-  if (!F2) return false;
   if ((theTria1->GetEntityType() == SMDSEntity_Triangle) &&
-      (theTria2->GetEntityType() == SMDSEntity_Triangle)) {
-
+      (theTria2->GetEntityType() == SMDSEntity_Triangle))
+  {
     //  1 +--+ A  theTria1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
     //    | /|    theTria2: ( B A 2 ) B->1 ( 1 A 2 )   |\ |
     //    |/ |                                         | \|
@@ -810,9 +805,9 @@ bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1,
     gp_Pnt xyz;
     if ( F.IsNull() )
     {
-      xyz = ( SMESH_TNodeXYZ( nodes[3] ) +
-              SMESH_TNodeXYZ( nodes[4] ) +
-              SMESH_TNodeXYZ( nodes[5] )) / 3.;
+      xyz = ( SMESH_NodeXYZ( nodes[3] ) +
+              SMESH_NodeXYZ( nodes[4] ) +
+              SMESH_NodeXYZ( nodes[5] )) / 3.;
     }
     else
     {
@@ -885,19 +880,16 @@ static bool findTriangles(const SMDS_MeshNode *    theNode1,
 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshNode * theNode1,
                                     const SMDS_MeshNode * theNode2)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
-
-  MESSAGE( "::InverseDiag()" );
+  ClearLastCreated();
 
   const SMDS_MeshElement *tr1, *tr2;
   if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
     return false;
 
-  const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( tr1 );
-  if (!F1) return false;
-  const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( tr2 );
-  if (!F2) return false;
+  if ( !dynamic_cast<const SMDS_MeshCell*>( tr1 ) ||
+       !dynamic_cast<const SMDS_MeshCell*>( tr2 ))
+    return false;
+
   if ((tr1->GetEntityType() == SMDSEntity_Triangle) &&
       (tr2->GetEntityType() == SMDSEntity_Triangle)) {
 
@@ -1008,37 +1000,33 @@ bool getQuadrangleNodes(const SMDS_MeshNode *    theQuadNodes [],
 bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
                                    const SMDS_MeshNode * theNode2)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
-
-  MESSAGE( "::DeleteDiag()" );
+  ClearLastCreated();
 
   const SMDS_MeshElement *tr1, *tr2;
   if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
     return false;
 
-  const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( tr1 );
-  if (!F1) return false;
-  const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( tr2 );
-  if (!F2) return false;
+  if ( !dynamic_cast<const SMDS_MeshCell*>( tr1 ) ||
+       !dynamic_cast<const SMDS_MeshCell*>( tr2 ))
+    return false;
+
   SMESHDS_Mesh * aMesh = GetMeshDS();
 
   if ((tr1->GetEntityType() == SMDSEntity_Triangle) &&
-      (tr2->GetEntityType() == SMDSEntity_Triangle)) {
-
+      (tr2->GetEntityType() == SMDSEntity_Triangle))
+  {
     const SMDS_MeshNode* aNodes [ 4 ];
     if ( ! getQuadrangleNodes( aNodes, theNode1, theNode2, tr1, tr2 ))
       return false;
 
     const SMDS_MeshElement* newElem = 0;
     newElem = aMesh->AddFace( aNodes[0], aNodes[1], aNodes[2], aNodes[3] );
-    myLastCreatedElems.Append(newElem);
+    myLastCreatedElems.push_back(newElem);
     AddToSameGroups( newElem, tr1, aMesh );
     int aShapeId = tr1->getshapeId();
     if ( aShapeId )
-      {
-        aMesh->SetMeshElementOnShape( newElem, aShapeId );
-      }
+      aMesh->SetMeshElementOnShape( newElem, aShapeId );
+
     aMesh->RemoveElement( tr1 );
     aMesh->RemoveElement( tr2 );
 
@@ -1082,13 +1070,13 @@ bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
   const SMDS_MeshElement* newElem = 0;
   newElem = aMesh->AddFace( aNodes[0], aNodes[1], aNodes[2], aNodes[3],
                             aNodes[4], aNodes[5], aNodes[6], aNodes[7]);
-  myLastCreatedElems.Append(newElem);
+  myLastCreatedElems.push_back(newElem);
   AddToSameGroups( newElem, tr1, aMesh );
   int aShapeId = tr1->getshapeId();
   if ( aShapeId )
-    {
-      aMesh->SetMeshElementOnShape( newElem, aShapeId );
-    }
+  {
+    aMesh->SetMeshElementOnShape( newElem, aShapeId );
+  }
   aMesh->RemoveElement( tr1 );
   aMesh->RemoveElement( tr2 );
 
@@ -1105,9 +1093,7 @@ bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
 
 bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
 {
-  MESSAGE("Reorient");
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
   if (!theElem)
     return false;
@@ -1122,8 +1108,7 @@ bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
   const SMDSAbs_EntityType geomType = theElem->GetEntityType();
   if ( geomType == SMDSEntity_Polyhedra ) // polyhedron
   {
-    const SMDS_VtkVolume* aPolyedre =
-      dynamic_cast<const SMDS_VtkVolume*>( theElem );
+    const SMDS_MeshVolume* aPolyedre = SMDS_Mesh::DownCast< SMDS_MeshVolume >( theElem );
     if (!aPolyedre) {
       MESSAGE("Warning: bad volumic element");
       return false;
@@ -1166,7 +1151,7 @@ bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
  * \brief Reorient faces.
  * \param theFaces - the faces to reorient. If empty the whole mesh is meant
  * \param theDirection - desired direction of normal of \a theFace
- * \param theFace - one of \a theFaces that sould be oriented according to
+ * \param theFace - one of \a theFaces that should be oriented according to
  *        \a theDirection and whose orientation defines orientation of other faces
  * \return number of reoriented faces.
  */
@@ -1181,7 +1166,7 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
 
   if ( theFaces.empty() )
   {
-    SMDS_FaceIteratorPtr fIt = GetMeshDS()->facesIterator(/*idInceasingOrder=*/true);
+    SMDS_FaceIteratorPtr fIt = GetMeshDS()->facesIterator(/*idInceasingOrder=true*/);
     while ( fIt->more() )
       theFaces.insert( theFaces.end(), fIt->next() );
   }
@@ -1276,7 +1261,7 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
         }
         if ( otherFace && otherFace != theFace)
         {
-          // link must be reverse in otherFace if orientation ot otherFace
+          // link must be reverse in otherFace if orientation to otherFace
           // is same as that of theFace
           if ( abs(nodeInd2-nodeInd1) == 1 ? nodeInd2 > nodeInd1 : nodeInd1 > nodeInd2 )
           {
@@ -1312,7 +1297,7 @@ int SMESH_MeshEditor::Reorient2DBy3D (TIDSortedElemSet & theFaces,
   if ( theFaces.empty() )
     faceIt = GetMeshDS()->elementsIterator( SMDSAbs_Face );
   else
-    faceIt = elemSetIterator( theFaces );
+    faceIt = SMESHUtils::elemSetIterator( theFaces );
 
   vector< const SMDS_MeshNode* > faceNodes;
   TIDSortedElemSet checkedVolumes;
@@ -1401,17 +1386,17 @@ static double getBadRate (const SMDS_MeshElement*               theElem,
 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet &                   theElems,
                                   SMESH::Controls::NumericalFunctorPtr theCrit)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
   if ( !theCrit.get() )
     return false;
 
-  SMESHDS_Mesh * aMesh = GetMeshDS();
-
+  SMESHDS_Mesh *       aMesh = GetMeshDS();
   Handle(Geom_Surface) surface;
   SMESH_MesherHelper   helper( *GetMesh() );
 
+  myLastCreatedElems.reserve( theElems.size() * 2 );
+
   TIDSortedElemSet::iterator itElem;
   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
   {
@@ -1438,7 +1423,7 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet &                   theElems,
     const SMDS_MeshElement* newElem1 = 0;
     const SMDS_MeshElement* newElem2 = 0;
 
-    if ( !elem->IsQuadratic() ) // split liner quadrangle
+    if ( !elem->IsQuadratic() ) // split linear quadrangle
     {
       // for MaxElementLength2D functor we return minimum diagonal for splitting,
       // because aBadRate1=2*len(diagonal 1-3); aBadRate2=2*len(diagonal 2-4)
@@ -1480,8 +1465,8 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet &                   theElems,
 
     // care of a new element
 
-    myLastCreatedElems.Append(newElem1);
-    myLastCreatedElems.Append(newElem2);
+    myLastCreatedElems.push_back(newElem1);
+    myLastCreatedElems.push_back(newElem2);
     AddToSameGroups( newElem1, elem, aMesh );
     AddToSameGroups( newElem2, elem, aMesh );
 
@@ -1498,27 +1483,27 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet &                   theElems,
 //=======================================================================
 /*!
  * \brief Split each of given quadrangles into 4 triangles.
- * \param theElems - The faces to be splitted. If empty all faces are split.
+ * \param theElems - The faces to be split. If empty all faces are split.
  */
 //=======================================================================
 
 void SMESH_MeshEditor::QuadTo4Tri (TIDSortedElemSet & theElems)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
+  myLastCreatedElems.reserve( theElems.size() * 4 );
 
   SMESH_MesherHelper helper( *GetMesh() );
   helper.SetElementsOnShape( true );
 
   SMDS_ElemIteratorPtr faceIt;
   if ( theElems.empty() ) faceIt = GetMeshDS()->elementsIterator(SMDSAbs_Face);
-  else                    faceIt = elemSetIterator( theElems );
+  else                    faceIt = SMESHUtils::elemSetIterator( theElems );
 
   bool   checkUV;
   gp_XY  uv [9]; uv[8] = gp_XY(0,0);
   gp_XYZ xyz[9];
   vector< const SMDS_MeshNode* > nodes;
-  SMESHDS_SubMesh*               subMeshDS;
+  SMESHDS_SubMesh*               subMeshDS = 0;
   TopoDS_Face                    F;
   Handle(Geom_Surface)           surface;
   TopLoc_Location                loc;
@@ -1569,7 +1554,7 @@ void SMESH_MeshEditor::QuadTo4Tri (TIDSortedElemSet & theElems)
       if ( F.IsNull() )
       {
         for ( ; iN < nodes.size(); ++iN )
-          xyz[ iN ] = SMESH_TNodeXYZ( nodes[ iN ] );
+          xyz[ iN ] = SMESH_NodeXYZ( nodes[ iN ] );
 
         for ( ; iN < 8; ++iN ) // mid-side points of a linear qudrangle
           xyz[ iN ] = 0.5 * ( xyz[ iN - 4 ] + xyz[( iN - 3 )%4 ] );
@@ -1596,25 +1581,25 @@ void SMESH_MeshEditor::QuadTo4Tri (TIDSortedElemSet & theElems)
 
       nCentral = helper.AddNode( xyz[8].X(), xyz[8].Y(), xyz[8].Z(), /*id=*/0,
                                  uv[8].X(), uv[8].Y() );
-      myLastCreatedNodes.Append( nCentral );
+      myLastCreatedNodes.push_back( nCentral );
     }
 
     // create 4 triangles
 
-    GetMeshDS()->RemoveFreeElement( quad, subMeshDS, /*fromGroups=*/false );
-    
     helper.SetIsQuadratic  ( nodes.size() > 4 );
     helper.SetIsBiQuadratic( nodes.size() == 9 );
     if ( helper.GetIsQuadratic() )
       helper.AddTLinks( static_cast< const SMDS_MeshFace*>( quad ));
 
+    GetMeshDS()->RemoveFreeElement( quad, subMeshDS, /*fromGroups=*/false );
+
     for ( int i = 0; i < 4; ++i )
     {
       SMDS_MeshElement* tria = helper.AddFace( nodes[ i ],
                                                nodes[(i+1)%4],
                                                nCentral );
       ReplaceElemInGroups( tria, quad, GetMeshDS() );
-      myLastCreatedElems.Append( tria );
+      myLastCreatedElems.push_back( tria );
     }
   }
 }
@@ -1627,8 +1612,7 @@ void SMESH_MeshEditor::QuadTo4Tri (TIDSortedElemSet & theElems)
 int SMESH_MeshEditor::BestSplit (const SMDS_MeshElement*              theQuad,
                                  SMESH::Controls::NumericalFunctorPtr theCrit)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
   if (!theCrit.get())
     return -1;
@@ -1912,6 +1896,7 @@ namespace
         break;
       case SMDSEntity_Penta:
       case SMDSEntity_Quad_Penta:
+      case SMDSEntity_BiQuad_Penta:
         connVariants = thePentaTo3; nbTet = 3; nbVariants = 6;
         break;
       default:
@@ -2130,7 +2115,7 @@ namespace
 
     // No adjacent prisms. Select a variant with a best aspect ratio.
 
-    double badness[2] = { 0, 0 };
+    double badness[2] = { 0., 0. };
     static SMESH::Controls::NumericalFunctorPtr aspectRatio( new SMESH::Controls::AspectRatio);
     const SMDS_MeshNode** nodes = vol.GetNodes();
     for ( int variant = 0; variant < nbVariants; ++variant )
@@ -2294,7 +2279,7 @@ void SMESH_MeshEditor::SplitVolumes (const TFacetOfElem & theElems,
       volTool.GetBaryCenter( bc[0], bc[1], bc[2] );
       SMDS_MeshNode* gcNode = helper.AddNode( bc[0], bc[1], bc[2] );
       nodes.push_back( gcNode );
-      newNodes.Append( gcNode );
+      newNodes.push_back( gcNode );
     }
     if ( !splitMethod._faceBaryNode.empty() )
     {
@@ -2308,7 +2293,7 @@ void SMESH_MeshEditor::SplitVolumes (const TFacetOfElem & theElems,
         if ( !f_n->second )
         {
           volTool.GetFaceBaryCenter( iF_n->first, bc[0], bc[1], bc[2] );
-          newNodes.Append( f_n->second = helper.AddNode( bc[0], bc[1], bc[2] ));
+          newNodes.push_back( f_n->second = helper.AddNode( bc[0], bc[1], bc[2] ));
         }
         nodes.push_back( iF_n->second = f_n->second );
       }
@@ -2319,18 +2304,18 @@ void SMESH_MeshEditor::SplitVolumes (const TFacetOfElem & theElems,
     const int* volConn = splitMethod._connectivity;
     if ( splitMethod._nbCorners == 4 ) // tetra
       for ( int i = 0; i < splitMethod._nbSplits; ++i, volConn += splitMethod._nbCorners )
-        newElems.Append( splitVols[ i ] = helper.AddVolume( nodes[ volConn[0] ],
-                                                            nodes[ volConn[1] ],
-                                                            nodes[ volConn[2] ],
-                                                            nodes[ volConn[3] ]));
+        newElems.push_back( splitVols[ i ] = helper.AddVolume( nodes[ volConn[0] ],
+                                                               nodes[ volConn[1] ],
+                                                               nodes[ volConn[2] ],
+                                                               nodes[ volConn[3] ]));
     else // prisms
       for ( int i = 0; i < splitMethod._nbSplits; ++i, volConn += splitMethod._nbCorners )
-        newElems.Append( splitVols[ i ] = helper.AddVolume( nodes[ volConn[0] ],
-                                                            nodes[ volConn[1] ],
-                                                            nodes[ volConn[2] ],
-                                                            nodes[ volConn[3] ],
-                                                            nodes[ volConn[4] ],
-                                                            nodes[ volConn[5] ]));
+        newElems.push_back( splitVols[ i ] = helper.AddVolume( nodes[ volConn[0] ],
+                                                               nodes[ volConn[1] ],
+                                                               nodes[ volConn[2] ],
+                                                               nodes[ volConn[3] ],
+                                                               nodes[ volConn[4] ],
+                                                               nodes[ volConn[5] ]));
 
     ReplaceElemInGroups( elem, splitVols, GetMeshDS() );
 
@@ -2374,7 +2359,7 @@ void SMESH_MeshEditor::SplitVolumes (const TFacetOfElem & theElems,
           if ( fSubMesh ) // update position of the bary node on geometry
           {
             if ( subMesh )
-              subMesh->RemoveNode( baryNode, false );
+              subMesh->RemoveNode( baryNode );
             GetMeshDS()->SetNodeOnFace( baryNode, fSubMesh->GetID() );
             const TopoDS_Shape& s = GetMeshDS()->IndexToShape( fSubMesh->GetID() );
             if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
@@ -2383,7 +2368,7 @@ void SMESH_MeshEditor::SplitVolumes (const TFacetOfElem & theElems,
               gp_XY uv( 1e100, 1e100 );
               double distXYZ[4];
               if ( !fHelper.CheckNodeUV( TopoDS::Face( s ), baryNode,
-                                        uv, /*tol=*/1e-7, /*force=*/true, distXYZ ) &&
+                                         uv, /*tol=*/1e-7, /*force=*/true, distXYZ ) &&
                    uv.X() < 1e100 )
               {
                 // node is too far from the surface
@@ -2396,7 +2381,7 @@ void SMESH_MeshEditor::SplitVolumes (const TFacetOfElem & theElems,
         }
         else
         {
-          // among possible triangles create ones discribed by split method
+          // among possible triangles create ones described by split method
           const int* nInd = volTool.GetFaceNodesIndices( iF );
           int nbVariants = ( nbNodes == 4 ? 2 : nbNodes );
           int iCom = 0; // common node of triangle faces to split into
@@ -2437,7 +2422,7 @@ void SMESH_MeshEditor::SplitVolumes (const TFacetOfElem & theElems,
           if ( !triangles[ i ]) continue;
           if ( fSubMesh )
             fSubMesh->AddElement( triangles[ i ]);
-          newElems.Append( triangles[ i ]);
+          newElems.push_back( triangles[ i ]);
         }
         ReplaceElemInGroups( face, triangles, GetMeshDS() );
         GetMeshDS()->RemoveFreeElement( face, fSubMesh, /*fromGroups=*/false );
@@ -2474,7 +2459,7 @@ void SMESH_MeshEditor::GetHexaFacetsToSplit( TIDSortedElemSet& theHexas,
                                              const gp_Ax1&     theFacetNormal,
                                              TFacetOfElem &    theFacets)
 {
-  #define THIS_METHOD "SMESH_MeshEditor::GetHexaFacetsToSplit(): "
+#define THIS_METHOD "SMESH_MeshEditor::GetHexaFacetsToSplit(): "
 
   // Find a hexa closest to the location of theFacetNormal
 
@@ -2522,7 +2507,7 @@ void SMESH_MeshEditor::GetHexaFacetsToSplit( TIDSortedElemSet& theHexas,
 
   // Fill theFacets starting from facetID of startHex
 
-  // facets used for seach of volumes adjacent to already treated ones
+  // facets used for searching of volumes adjacent to already treated ones
   typedef pair< TFacetOfElem::iterator, int > TElemFacets;
   typedef map< TVolumeFaceKey, TElemFacets  > TFacetMap;
   TFacetMap facetsToCheck;
@@ -2928,13 +2913,10 @@ void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement*
 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
                                   const bool         the13Diag)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
-
-  MESSAGE( "::QuadToTri()" );
-
-  SMESHDS_Mesh * aMesh = GetMeshDS();
+  ClearLastCreated();
+  myLastCreatedElems.reserve( theElems.size() * 2 );
 
+  SMESHDS_Mesh *       aMesh = GetMeshDS();
   Handle(Geom_Surface) surface;
   SMESH_MesherHelper   helper( *GetMesh() );
 
@@ -2964,8 +2946,8 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
         newElem1 = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
         newElem2 = aMesh->AddFace( aNodes[3], aNodes[1], aNodes[2] );
       }
-      myLastCreatedElems.Append(newElem1);
-      myLastCreatedElems.Append(newElem2);
+      myLastCreatedElems.push_back(newElem1);
+      myLastCreatedElems.push_back(newElem2);
       // put a new triangle on the same shape and add to the same groups
       if ( aShapeId )
       {
@@ -2979,7 +2961,7 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
 
     // Quadratic quadrangle
 
-    else if ( elem->NbNodes() == 8 )
+    else if ( elem->NbNodes() >= 8 )
     {
       // get surface elem is on
       int aShapeId = FindShape( elem );
@@ -2996,55 +2978,37 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
         }
       }
 
-      const SMDS_MeshNode* aNodes [8];
-      const SMDS_MeshNode* inFaceNode = 0;
+      const SMDS_MeshNode* aNodes [9]; aNodes[8] = 0;
       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
-      int i = 0;
-      if ( helper.GetNodeUVneedInFaceNode() )
-        while ( itN->more() && !inFaceNode ) {
-          aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
-          if ( aNodes[ i-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
-          {
-            inFaceNode = aNodes[ i-1 ];
-          }
-        }
+      for ( int i = 0; itN->more(); ++i )
+        aNodes[ i ] = static_cast<const SMDS_MeshNode*>( itN->next() );
 
-      // find middle point for (0,1,2,3)
-      // and create a node in this point;
-      gp_XYZ p( 0,0,0 );
-      if ( surface.IsNull() ) {
-        for(i=0; i<4; i++)
-          p += gp_XYZ(aNodes[i]->X(), aNodes[i]->Y(), aNodes[i]->Z() );
-        p /= 4;
-      }
-      else {
-        TopoDS_Face geomFace = TopoDS::Face( helper.GetSubShape() );
-        gp_XY uv( 0,0 );
-        for(i=0; i<4; i++)
-          uv += helper.GetNodeUV( geomFace, aNodes[i], inFaceNode );
-        uv /= 4.;
-        p = surface->Value( uv.X(), uv.Y() ).XYZ();
+      const SMDS_MeshNode* centrNode = aNodes[8];
+      if ( centrNode == 0 )
+      {
+        centrNode = helper.GetCentralNode( aNodes[0], aNodes[1], aNodes[2], aNodes[3],
+                                           aNodes[4], aNodes[5], aNodes[6], aNodes[7],
+                                           surface.IsNull() );
+        myLastCreatedNodes.push_back(centrNode);
       }
-      const SMDS_MeshNode* newN = aMesh->AddNode( p.X(), p.Y(), p.Z() );
-      myLastCreatedNodes.Append(newN);
 
       // create a new element
       const SMDS_MeshElement* newElem1 = 0;
       const SMDS_MeshElement* newElem2 = 0;
       if ( the13Diag ) {
         newElem1 = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
-                                  aNodes[6], aNodes[7], newN );
+                                  aNodes[6], aNodes[7], centrNode );
         newElem2 = aMesh->AddFace(aNodes[2], aNodes[0], aNodes[1],
-                                  newN,      aNodes[4], aNodes[5] );
+                                  centrNode, aNodes[4], aNodes[5] );
       }
       else {
         newElem1 = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
-                                  aNodes[7], aNodes[4], newN );
+                                  aNodes[7], aNodes[4], centrNode );
         newElem2 = aMesh->AddFace(aNodes[3], aNodes[1], aNodes[2],
-                                  newN,      aNodes[5], aNodes[6] );
+                                  centrNode, aNodes[5], aNodes[6] );
       }
-      myLastCreatedElems.Append(newElem1);
-      myLastCreatedElems.Append(newElem2);
+      myLastCreatedElems.push_back(newElem1);
+      myLastCreatedElems.push_back(newElem2);
       // put a new triangle on the same shape and add to the same groups
       if ( aShapeId )
       {
@@ -3168,10 +3132,8 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
                                   SMESH::Controls::NumericalFunctorPtr theCrit,
                                   const double                         theMaxAngle)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
-
-  MESSAGE( "::TriToQuad()" );
+  ClearLastCreated();
+  myLastCreatedElems.reserve( theElems.size() / 2 );
 
   if ( !theCrit.get() )
     return false;
@@ -3277,7 +3239,7 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
       if ( startElem ) {
         // Get candidates to be fused
         const SMDS_MeshElement *tr1 = startElem, *tr2 = 0, *tr3 = 0;
-        const SMESH_TLink *link12, *link13;
+        const SMESH_TLink *link12 = 0, *link13 = 0;
         startElem = 0;
         ASSERT( mapEl_setLi.find( tr1 ) != mapEl_setLi.end() );
         set< SMESH_TLink >& setLi = mapEl_setLi[ tr1 ];
@@ -3355,7 +3317,7 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
           {
             const SMDS_MeshElement* newElem = 0;
             newElem = aMesh->AddFace(n12[0], n12[1], n12[2], n12[3] );
-            myLastCreatedElems.Append(newElem);
+            myLastCreatedElems.push_back(newElem);
             AddToSameGroups( newElem, tr1, aMesh );
             int aShapeId = tr1->getshapeId();
             if ( aShapeId )
@@ -3386,7 +3348,7 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
             else
               newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3],
                                        aNodes[4], aNodes[5], aNodes[6], aNodes[7]);
-            myLastCreatedElems.Append(newElem);
+            myLastCreatedElems.push_back(newElem);
             AddToSameGroups( newElem, tr1, aMesh );
             int aShapeId = tr1->getshapeId();
             if ( aShapeId )
@@ -3409,7 +3371,7 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
           if ( tr1->NbNodes() == 3 ) {
             const SMDS_MeshElement* newElem = 0;
             newElem = aMesh->AddFace(n13[0], n13[1], n13[2], n13[3] );
-            myLastCreatedElems.Append(newElem);
+            myLastCreatedElems.push_back(newElem);
             AddToSameGroups( newElem, tr1, aMesh );
             int aShapeId = tr1->getshapeId();
             if ( aShapeId )
@@ -3440,7 +3402,7 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
             else
               newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3],
                                        aNodes[4], aNodes[5], aNodes[6], aNodes[7]);
-            myLastCreatedElems.Append(newElem);
+            myLastCreatedElems.push_back(newElem);
             AddToSameGroups( newElem, tr1, aMesh );
             int aShapeId = tr1->getshapeId();
             if ( aShapeId )
@@ -3468,286 +3430,6 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
   return true;
 }
 
-
-/*#define DUMPSO(txt) \
-//  cout << txt << endl;
-//=============================================================================
-//
-//
-//
-//=============================================================================
-static void swap( int i1, int i2, int idNodes[], gp_Pnt P[] )
-{
-if ( i1 == i2 )
-return;
-int tmp = idNodes[ i1 ];
-idNodes[ i1 ] = idNodes[ i2 ];
-idNodes[ i2 ] = tmp;
-gp_Pnt Ptmp = P[ i1 ];
-P[ i1 ] = P[ i2 ];
-P[ i2 ] = Ptmp;
-DUMPSO( i1 << "(" << idNodes[ i2 ] << ") <-> " << i2 << "(" << idNodes[ i1 ] << ")");
-}
-
-//=======================================================================
-//function : SortQuadNodes
-//purpose  : Set 4 nodes of a quadrangle face in a good order.
-//           Swap 1<->2 or 2<->3 nodes and correspondingly return
-//           1 or 2 else 0.
-//=======================================================================
-
-int SMESH_MeshEditor::SortQuadNodes (const SMDS_Mesh * theMesh,
-int               idNodes[] )
-{
-  gp_Pnt P[4];
-  int i;
-  for ( i = 0; i < 4; i++ ) {
-    const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
-    if ( !n ) return 0;
-    P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
-  }
-
-  gp_Vec V1(P[0], P[1]);
-  gp_Vec V2(P[0], P[2]);
-  gp_Vec V3(P[0], P[3]);
-
-  gp_Vec Cross1 = V1 ^ V2;
-  gp_Vec Cross2 = V2 ^ V3;
-
-  i = 0;
-  if (Cross1.Dot(Cross2) < 0)
-  {
-    Cross1 = V2 ^ V1;
-    Cross2 = V1 ^ V3;
-
-    if (Cross1.Dot(Cross2) < 0)
-      i = 2;
-    else
-      i = 1;
-    swap ( i, i + 1, idNodes, P );
-
-    //     for ( int ii = 0; ii < 4; ii++ ) {
-    //       const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
-    //       DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
-    //     }
-  }
-  return i;
-}
-
-//=======================================================================
-//function : SortHexaNodes
-//purpose  : Set 8 nodes of a hexahedron in a good order.
-//           Return success status
-//=======================================================================
-
-bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh,
-                                      int               idNodes[] )
-{
-  gp_Pnt P[8];
-  int i;
-  DUMPSO( "INPUT: ========================================");
-  for ( i = 0; i < 8; i++ ) {
-    const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
-    if ( !n ) return false;
-    P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
-    DUMPSO( i << "(" << idNodes[i] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
-  }
-  DUMPSO( "========================================");
-
-
-  set<int> faceNodes;  // ids of bottom face nodes, to be found
-  set<int> checkedId1; // ids of tried 2-nd nodes
-  Standard_Real leastDist = DBL_MAX; // dist of the 4-th node from 123 plane
-  const Standard_Real tol = 1.e-6;   // tolerance to find nodes in plane
-  int iMin, iLoop1 = 0;
-
-  // Loop to try the 2-nd nodes
-
-  while ( leastDist > DBL_MIN && ++iLoop1 < 8 )
-  {
-    // Find not checked 2-nd node
-    for ( i = 1; i < 8; i++ )
-      if ( checkedId1.find( idNodes[i] ) == checkedId1.end() ) {
-        int id1 = idNodes[i];
-        swap ( 1, i, idNodes, P );
-        checkedId1.insert ( id1 );
-        break;
-      }
-
-    // Find the 3-d node so that 1-2-3 triangle to be on a hexa face,
-    // ie that all but meybe one (id3 which is on the same face) nodes
-    // lay on the same side from the triangle plane.
-
-    bool manyInPlane = false; // more than 4 nodes lay in plane
-    int iLoop2 = 0;
-    while ( ++iLoop2 < 6 ) {
-
-      // get 1-2-3 plane coeffs
-      Standard_Real A, B, C, D;
-      gp_Vec N = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
-      if ( N.SquareMagnitude() > gp::Resolution() )
-      {
-        gp_Pln pln ( P[0], N );
-        pln.Coefficients( A, B, C, D );
-
-        // find the node (iMin) closest to pln
-        Standard_Real dist[ 8 ], minDist = DBL_MAX;
-        set<int> idInPln;
-        for ( i = 3; i < 8; i++ ) {
-          dist[i] = A * P[i].X() + B * P[i].Y() + C * P[i].Z() + D;
-          if ( fabs( dist[i] ) < minDist ) {
-            minDist = fabs( dist[i] );
-            iMin = i;
-          }
-          if ( fabs( dist[i] ) <= tol )
-            idInPln.insert( idNodes[i] );
-        }
-
-        // there should not be more than 4 nodes in bottom plane
-        if ( idInPln.size() > 1 )
-        {
-          DUMPSO( "### idInPln.size() = " << idInPln.size());
-          // idInPlane does not contain the first 3 nodes
-          if ( manyInPlane || idInPln.size() == 5)
-            return false; // all nodes in one plane
-          manyInPlane = true;
-
-          // set the 1-st node to be not in plane
-          for ( i = 3; i < 8; i++ ) {
-            if ( idInPln.find( idNodes[ i ] ) == idInPln.end() ) {
-              DUMPSO( "### Reset 0-th node");
-              swap( 0, i, idNodes, P );
-              break;
-            }
-          }
-
-          // reset to re-check second nodes
-          leastDist = DBL_MAX;
-          faceNodes.clear();
-          checkedId1.clear();
-          iLoop1 = 0;
-          break; // from iLoop2;
-        }
-
-        // check that the other 4 nodes are on the same side
-        bool sameSide = true;
-        bool isNeg = dist[ iMin == 3 ? 4 : 3 ] <= 0.;
-        for ( i = 3; sameSide && i < 8; i++ ) {
-          if ( i != iMin )
-            sameSide = ( isNeg == dist[i] <= 0.);
-        }
-
-        // keep best solution
-        if ( sameSide && minDist < leastDist ) {
-          leastDist = minDist;
-          faceNodes.clear();
-          faceNodes.insert( idNodes[ 1 ] );
-          faceNodes.insert( idNodes[ 2 ] );
-          faceNodes.insert( idNodes[ iMin ] );
-          DUMPSO( "loop " << iLoop2 << " id2 " << idNodes[ 1 ] << " id3 " << idNodes[ 2 ]
-                  << " leastDist = " << leastDist);
-          if ( leastDist <= DBL_MIN )
-            break;
-        }
-      }
-
-      // set next 3-d node to check
-      int iNext = 2 + iLoop2;
-      if ( iNext < 8 ) {
-        DUMPSO( "Try 2-nd");
-        swap ( 2, iNext, idNodes, P );
-      }
-    } // while ( iLoop2 < 6 )
-  } // iLoop1
-
-  if ( faceNodes.empty() ) return false;
-
-  // Put the faceNodes in proper places
-  for ( i = 4; i < 8; i++ ) {
-    if ( faceNodes.find( idNodes[ i ] ) != faceNodes.end() ) {
-      // find a place to put
-      int iTo = 1;
-      while ( faceNodes.find( idNodes[ iTo ] ) != faceNodes.end() )
-        iTo++;
-      DUMPSO( "Set faceNodes");
-      swap ( iTo, i, idNodes, P );
-    }
-  }
-
-
-  // Set nodes of the found bottom face in good order
-  DUMPSO( " Found bottom face: ");
-  i = SortQuadNodes( theMesh, idNodes );
-  if ( i ) {
-    gp_Pnt Ptmp = P[ i ];
-    P[ i ] = P[ i+1 ];
-    P[ i+1 ] = Ptmp;
-  }
-  //   else
-  //     for ( int ii = 0; ii < 4; ii++ ) {
-  //       const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
-  //       DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
-  //    }
-
-  // Gravity center of the top and bottom faces
-  gp_Pnt aGCb = ( P[0].XYZ() + P[1].XYZ() + P[2].XYZ() + P[3].XYZ() ) / 4.;
-  gp_Pnt aGCt = ( P[4].XYZ() + P[5].XYZ() + P[6].XYZ() + P[7].XYZ() ) / 4.;
-
-  // Get direction from the bottom to the top face
-  gp_Vec upDir ( aGCb, aGCt );
-  Standard_Real upDirSize = upDir.Magnitude();
-  if ( upDirSize <= gp::Resolution() ) return false;
-  upDir / upDirSize;
-
-  // Assure that the bottom face normal points up
-  gp_Vec Nb = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
-  Nb += gp_Vec (P[0], P[2]).Crossed( gp_Vec (P[0], P[3]) );
-  if ( Nb.Dot( upDir ) < 0 ) {
-    DUMPSO( "Reverse bottom face");
-    swap( 1, 3, idNodes, P );
-  }
-
-  // Find 5-th node - the one closest to the 1-st among the last 4 nodes.
-  Standard_Real minDist = DBL_MAX;
-  for ( i = 4; i < 8; i++ ) {
-    // projection of P[i] to the plane defined by P[0] and upDir
-    gp_Pnt Pp = P[i].Translated( upDir * ( upDir.Dot( gp_Vec( P[i], P[0] ))));
-    Standard_Real sqDist = P[0].SquareDistance( Pp );
-    if ( sqDist < minDist ) {
-      minDist = sqDist;
-      iMin = i;
-    }
-  }
-  DUMPSO( "Set 4-th");
-  swap ( 4, iMin, idNodes, P );
-
-  // Set nodes of the top face in good order
-  DUMPSO( "Sort top face");
-  i = SortQuadNodes( theMesh, &idNodes[4] );
-  if ( i ) {
-    i += 4;
-    gp_Pnt Ptmp = P[ i ];
-    P[ i ] = P[ i+1 ];
-    P[ i+1 ] = Ptmp;
-  }
-
-  // Assure that direction of the top face normal is from the bottom face
-  gp_Vec Nt = gp_Vec (P[4], P[5]).Crossed( gp_Vec (P[4], P[6]) );
-  Nt += gp_Vec (P[4], P[6]).Crossed( gp_Vec (P[4], P[7]) );
-  if ( Nt.Dot( upDir ) < 0 ) {
-    DUMPSO( "Reverse top face");
-    swap( 5, 7, idNodes, P );
-  }
-
-  //   DUMPSO( "OUTPUT: ========================================");
-  //   for ( i = 0; i < 8; i++ ) {
-  //     float *p = ugrid->GetPoint(idNodes[i]);
-  //     DUMPSO( i << "(" << idNodes[i] << ") : " << p[0] << " " << p[1] << " " << p[2]);
-  //   }
-
-  return true;
-}*/
-
 //================================================================================
 /*!
  * \brief Return nodes linked to the given one
@@ -3953,10 +3635,7 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
                                double                      theTgtAspectRatio,
                                const bool                  the2D)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
-
-  MESSAGE((theSmoothMethod==LAPLACIAN ? "LAPLACIAN" : "CENTROIDAL") << "--::Smooth()");
+  ClearLastCreated();
 
   if ( theTgtAspectRatio < 1.0 )
     theTgtAspectRatio = 1.0;
@@ -4016,11 +3695,11 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
       fToler2 = BRep_Tool::Tolerance( face );
       fToler2 *= fToler2 * 10.;
       isUPeriodic = surface->IsUPeriodic();
-      if ( isUPeriodic )
-        surface->UPeriod();
+      // if ( isUPeriodic )
+      //   surface->UPeriod();
       isVPeriodic = surface->IsVPeriodic();
-      if ( isVPeriodic )
-        surface->VPeriod();
+      // if ( isVPeriodic )
+      //   surface->VPeriod();
       surface->Bounds( u1, u2, v1, v2 );
       helper.SetSubShape( face );
     }
@@ -4079,9 +3758,9 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
         {
           // check if all faces around the node are on faceSubMesh
           // because a node on edge may be bound to face
-          SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Face);
           bool all = true;
           if ( faceSubMesh ) {
+            SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Face);
             while ( eIt->more() && all ) {
               const SMDS_MeshElement* e = eIt->next();
               all = faceSubMesh->Contains( e );
@@ -4149,7 +3828,7 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
         // }
         if ( project ) { // compute new UV
           gp_XY newUV;
-          gp_Pnt pNode = SMESH_TNodeXYZ( node );
+          gp_Pnt pNode = SMESH_NodeXYZ( node );
           if ( !getClosestUV( projector, pNode, newUV )) {
             MESSAGE("Node Projection Failed " << node);
           }
@@ -4162,7 +3841,7 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
             // if ( posType != SMDS_TOP_3DSPACE )
             //   dist2 = pNode.SquareDistance( surface->Value( newUV.X(), newUV.Y() ));
             // if ( dist2 < dist1 )
-              uv = newUV;
+            uv = newUV;
           }
         }
         // store UV in the map
@@ -4364,18 +4043,18 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
         }
       }
       if ( maxRatio <= theTgtAspectRatio ) {
-        MESSAGE("-- quality achived --");
+        //MESSAGE("-- quality achieved --");
         break;
       }
       if (it+1 == theNbIterations) {
-        MESSAGE("-- Iteration limit exceeded --");
+        //MESSAGE("-- Iteration limit exceeded --");
       }
     } // smoothing iterations
 
-    MESSAGE(" Face id: " << *fId <<
-            " Nb iterstions: " << it <<
-            " Displacement: " << maxDisplacement <<
-            " Aspect Ratio " << maxRatio);
+    // MESSAGE(" Face id: " << *fId <<
+    //         " Nb iterstions: " << it <<
+    //         " Displacement: " << maxDisplacement <<
+    //         " Aspect Ratio " << maxRatio);
 
     // ---------------------------------------
     // new nodes positions are computed,
@@ -4404,7 +4083,7 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
         const SMDS_MeshElement* QF = *elemIt;
         if ( QF->IsQuadratic() )
         {
-          nodes.assign( SMDS_MeshElement::iterator( QF->interlacedNodesElemIterator() ),
+          nodes.assign( SMDS_MeshElement::iterator( QF->interlacedNodesIterator() ),
                         SMDS_MeshElement::iterator() );
           nodes.push_back( nodes[0] );
           gp_Pnt xyz;
@@ -4418,9 +4097,9 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
               xyz = surface->Value( uv.X(), uv.Y() );
             }
             else {
-              xyz = 0.5 * ( SMESH_TNodeXYZ( nodes[i-1] ) + SMESH_TNodeXYZ( nodes[i+1] ));
+              xyz = 0.5 * ( SMESH_NodeXYZ( nodes[i-1] ) + SMESH_NodeXYZ( nodes[i+1] ));
             }
-            if (( SMESH_TNodeXYZ( nodes[i] ) - xyz.XYZ() ).Modulus() > disttol )
+            if (( SMESH_NodeXYZ( nodes[i] ) - xyz.XYZ() ).Modulus() > disttol )
               // we have to move a medium node
               aMesh->MoveNode( nodes[i], xyz.X(), xyz.Y(), xyz.Z() );
           }
@@ -4448,8 +4127,8 @@ namespace
                  const int                           iNotSame)
   {
 
-    SMESH_TNodeXYZ pP = prevNodes[ iNotSame ];
-    SMESH_TNodeXYZ pN = nextNodes[ iNotSame ];
+    SMESH_NodeXYZ pP = prevNodes[ iNotSame ];
+    SMESH_NodeXYZ pN = nextNodes[ iNotSame ];
     gp_XYZ extrDir( pN - pP ), faceNorm;
     SMESH_MeshAlgos::FaceNormal( face, faceNorm, /*normalized=*/false );
 
@@ -4494,7 +4173,6 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
                                     const size_t                          nbSteps,
                                     SMESH_SequenceOfElemPtr&              srcElements)
 {
-  //MESSAGE("sweepElement " << nbSteps);
   SMESHDS_Mesh* aMesh = GetMeshDS();
 
   const int           nbNodes = elem->NbNodes();
@@ -4591,7 +4269,7 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
         std::swap( itNN[0],    itNN[1] );
         std::swap( prevNod[0], prevNod[1] );
         std::swap( nextNod[0], nextNod[1] );
-        isSingleNode.swap( isSingleNode[0], isSingleNode[1] );
+        std::swap( isSingleNode[0], isSingleNode[1] );
         if ( nbSame > 0 )
           sames[0] = 1 - sames[0];
         iNotSameNode = 1 - iNotSameNode;
@@ -4683,55 +4361,55 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
         break;
       }
       case SMDSEntity_Triangle: // TRIANGLE --->
-        {
-          if ( nbDouble > 0 ) break;
-          if ( nbSame == 0 )       // ---> pentahedron
-            aNewElem = aMesh->AddVolume (prevNod[ 0 ], prevNod[ 1 ], prevNod[ 2 ],
-                                         nextNod[ 0 ], nextNod[ 1 ], nextNod[ 2 ] );
-
-          else if ( nbSame == 1 )  // ---> pyramid
-            aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iAfterSame ],
-                                         nextNod[ iAfterSame ],  nextNod[ iBeforeSame ],
-                                         nextNod[ iSameNode ]);
-
-          else // 2 same nodes:       ---> tetrahedron
-            aNewElem = aMesh->AddVolume (prevNod[ 0 ], prevNod[ 1 ], prevNod[ 2 ],
-                                         nextNod[ iNotSameNode ]);
-          break;
-        }
+      {
+        if ( nbDouble > 0 ) break;
+        if ( nbSame == 0 )       // ---> pentahedron
+          aNewElem = aMesh->AddVolume (prevNod[ 0 ], prevNod[ 1 ], prevNod[ 2 ],
+                                       nextNod[ 0 ], nextNod[ 1 ], nextNod[ 2 ] );
+
+        else if ( nbSame == 1 )  // ---> pyramid
+          aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iAfterSame ],
+                                       nextNod[ iAfterSame ],  nextNod[ iBeforeSame ],
+                                       nextNod[ iSameNode ]);
+
+        else // 2 same nodes:       ---> tetrahedron
+          aNewElem = aMesh->AddVolume (prevNod[ 0 ], prevNod[ 1 ], prevNod[ 2 ],
+                                       nextNod[ iNotSameNode ]);
+        break;
+      }
       case SMDSEntity_Quad_Edge: // sweep quadratic EDGE --->
+      {
+        if ( nbSame == 2 )
+          return;
+        if ( nbDouble+nbSame == 2 )
         {
-          if ( nbSame == 2 )
-            return;
-          if ( nbDouble+nbSame == 2 )
-          {
-            if(nbSame==0) {      // ---> quadratic quadrangle
-              aNewElem = aMesh->AddFace(prevNod[0], prevNod[1], nextNod[1], nextNod[0],
-                                        prevNod[2], midlNod[1], nextNod[2], midlNod[0]);
-            }
-            else { //(nbSame==1) // ---> quadratic triangle
-              if(sames[0]==2) {
-                return; // medium node on axis
-              }
-              else if(sames[0]==0)
-                aNewElem = aMesh->AddFace(prevNod[0], prevNod[1], nextNod[1],
-                                          prevNod[2], midlNod[1], nextNod[2] );
-              else // sames[0]==1
-                aNewElem = aMesh->AddFace(prevNod[0], prevNod[1], nextNod[0],
-                                          prevNod[2], nextNod[2], midlNod[0]);
-            }
+          if(nbSame==0) {      // ---> quadratic quadrangle
+            aNewElem = aMesh->AddFace(prevNod[0], prevNod[1], nextNod[1], nextNod[0],
+                                      prevNod[2], midlNod[1], nextNod[2], midlNod[0]);
           }
-          else if ( nbDouble == 3 )
-          {
-            if ( nbSame == 0 ) {  // ---> bi-quadratic quadrangle
-              aNewElem = aMesh->AddFace(prevNod[0], prevNod[1], nextNod[1], nextNod[0],
-                                        prevNod[2], midlNod[1], nextNod[2], midlNod[0], midlNod[2]);
+          else { //(nbSame==1) // ---> quadratic triangle
+            if(sames[0]==2) {
+              return; // medium node on axis
             }
+            else if(sames[0]==0)
+              aNewElem = aMesh->AddFace(prevNod[0], prevNod[1], nextNod[1],
+                                        prevNod[2], midlNod[1], nextNod[2] );
+            else // sames[0]==1
+              aNewElem = aMesh->AddFace(prevNod[0], prevNod[1], nextNod[0],
+                                        prevNod[2], nextNod[2], midlNod[0]);
           }
-          else
-            return;
-          break;
         }
+        else if ( nbDouble == 3 )
+        {
+          if ( nbSame == 0 ) {  // ---> bi-quadratic quadrangle
+            aNewElem = aMesh->AddFace(prevNod[0], prevNod[1], nextNod[1], nextNod[0],
+                                      prevNod[2], midlNod[1], nextNod[2], midlNod[0], midlNod[2]);
+          }
+        }
+        else
+          return;
+        break;
+      }
       case SMDSEntity_Quadrangle: { // sweep QUADRANGLE --->
         if ( nbDouble > 0 ) break;
 
@@ -4953,8 +4631,8 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
 
     if ( aNewElem ) {
       newElems.push_back( aNewElem );
-      myLastCreatedElems.Append(aNewElem);
-      srcElements.Append( elem );
+      myLastCreatedElems.push_back(aNewElem);
+      srcElements.push_back( elem );
     }
 
     // set new prev nodes
@@ -5044,19 +4722,19 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
       if ( !isQuadratic ) {
         if ( !aMesh->FindEdge( vecNewNodes[ 0 ]->second.back(),
                                vecNewNodes[ 1 ]->second.back())) {
-          myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
-                                                   vecNewNodes[ 1 ]->second.back()));
-          srcElements.Append( elem );
+          myLastCreatedElems.push_back(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
+                                                      vecNewNodes[ 1 ]->second.back()));
+          srcElements.push_back( elem );
         }
       }
       else {
         if ( !aMesh->FindEdge( vecNewNodes[ 0 ]->second.back(),
                                vecNewNodes[ 1 ]->second.back(),
                                vecNewNodes[ 2 ]->second.back())) {
-          myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
-                                                   vecNewNodes[ 1 ]->second.back(),
-                                                   vecNewNodes[ 2 ]->second.back()));
-          srcElements.Append( elem );
+          myLastCreatedElems.push_back(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
+                                                      vecNewNodes[ 1 ]->second.back(),
+                                                      vecNewNodes[ 2 ]->second.back()));
+          srcElements.push_back( elem );
         }
       }
     }
@@ -5084,14 +4762,14 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
           // make a new edge and a ceiling for a new edge
           const SMDS_MeshElement* edge;
           if ( ! ( edge = aMesh->FindEdge( n1, n2 ))) {
-            myLastCreatedElems.Append( edge = aMesh->AddEdge( n1, n2 )); // free link edge
-            srcElements.Append( myLastCreatedElems.Last() );
+            myLastCreatedElems.push_back( edge = aMesh->AddEdge( n1, n2 )); // free link edge
+            srcElements.push_back( myLastCreatedElems.back() );
           }
           n1 = vecNewNodes[ iNode ]->second.back();
           n2 = vecNewNodes[ iNext ]->second.back();
           if ( !aMesh->FindEdge( n1, n2 )) {
-            myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // new edge ceiling
-            srcElements.Append( edge );
+            myLastCreatedElems.push_back(aMesh->AddEdge( n1, n2 )); // new edge ceiling
+            srcElements.push_back( edge );
           }
         }
       }
@@ -5112,15 +4790,15 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
           // make an edge and a ceiling for a new edge
           // find medium node
           if ( !aMesh->FindEdge( n1, n2, n3 )) {
-            myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // free link edge
-            srcElements.Append( elem );
+            myLastCreatedElems.push_back(aMesh->AddEdge( n1, n2, n3 )); // free link edge
+            srcElements.push_back( elem );
           }
           n1 = vecNewNodes[ iNode ]->second.back();
           n2 = vecNewNodes[ iNext ]->second.back();
           n3 = vecNewNodes[ iNode+nbn ]->second.back();
           if ( !aMesh->FindEdge( n1, n2, n3 )) {
-            myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // ceiling edge
-            srcElements.Append( elem );
+            myLastCreatedElems.push_back(aMesh->AddEdge( n1, n2, n3 )); // ceiling edge
+            srcElements.push_back( elem );
           }
         }
       }
@@ -5154,9 +4832,9 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
         SMDS_VolumeTool vTool( *v, /*ignoreCentralNodes=*/false );
         int iF, nbF = vTool.NbFaces();
         for ( iF = 0; iF < nbF; iF ++ ) {
-          if (vTool.IsFreeFace( iF ) &&
-              vTool.GetFaceNodes( iF, faceNodeSet ) &&
-              initNodeSet != faceNodeSet) // except an initial face
+          if ( vTool.IsFreeFace( iF ) &&
+               vTool.GetFaceNodes( iF, faceNodeSet ) &&
+               initNodeSet != faceNodeSet) // except an initial face
           {
             if ( nbSteps == 1 && faceNodeSet == topNodeSet )
               continue;
@@ -5213,8 +4891,8 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
                 if ( f )
                   aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
                 else
-                  myLastCreatedElems.Append(aMesh->AddFace( newOrder[ 0 ], newOrder[ 1 ],
-                                                            newOrder[ 2 ] ));
+                  myLastCreatedElems.push_back(aMesh->AddFace( newOrder[ 0 ], newOrder[ 1 ],
+                                                               newOrder[ 2 ] ));
               }
             }
             else if ( nbn == 4 )       ///// quadrangle
@@ -5228,8 +4906,8 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
                 if ( f )
                   aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
                 else
-                  myLastCreatedElems.Append(aMesh->AddFace( newOrder[ 0 ], newOrder[ 1 ],
-                                                            newOrder[ 2 ], newOrder[ 3 ]));
+                  myLastCreatedElems.push_back(aMesh->AddFace( newOrder[ 0 ], newOrder[ 1 ],
+                                                               newOrder[ 2 ], newOrder[ 3 ]));
               }
             }
             else if ( nbn == 6 && isQuadratic ) /////// quadratic triangle
@@ -5247,12 +4925,12 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
                 if ( f )
                   aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
                 else
-                  myLastCreatedElems.Append(aMesh->AddFace( newOrder[ 0 ],
-                                                            newOrder[ 1 ],
-                                                            newOrder[ 2 ],
-                                                            newOrder[ 3 ],
-                                                            newOrder[ 4 ],
-                                                            newOrder[ 5 ] ));
+                  myLastCreatedElems.push_back(aMesh->AddFace( newOrder[ 0 ],
+                                                               newOrder[ 1 ],
+                                                               newOrder[ 2 ],
+                                                               newOrder[ 3 ],
+                                                               newOrder[ 4 ],
+                                                               newOrder[ 5 ] ));
               }
             }
             else if ( nbn == 8 && isQuadratic ) /////// quadratic quadrangle
@@ -5273,10 +4951,10 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
                 if ( f )
                   aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
                 else
-                  myLastCreatedElems.Append(aMesh->AddFace(newOrder[ 0 ], newOrder[ 1 ],
-                                                           newOrder[ 2 ], newOrder[ 3 ],
-                                                           newOrder[ 4 ], newOrder[ 5 ],
-                                                           newOrder[ 6 ], newOrder[ 7 ]));
+                  myLastCreatedElems.push_back(aMesh->AddFace(newOrder[ 0 ], newOrder[ 1 ],
+                                                              newOrder[ 2 ], newOrder[ 3 ],
+                                                              newOrder[ 4 ], newOrder[ 5 ],
+                                                              newOrder[ 6 ], newOrder[ 7 ]));
               }
             }
             else if ( nbn == 9 && isQuadratic ) /////// bi-quadratic quadrangle
@@ -5298,11 +4976,11 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
                 if ( f )
                   aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
                 else
-                  myLastCreatedElems.Append(aMesh->AddFace(newOrder[ 0 ], newOrder[ 1 ],
-                                                           newOrder[ 2 ], newOrder[ 3 ],
-                                                           newOrder[ 4 ], newOrder[ 5 ],
-                                                           newOrder[ 6 ], newOrder[ 7 ],
-                                                           newOrder[ 8 ]));
+                  myLastCreatedElems.push_back(aMesh->AddFace(newOrder[ 0 ], newOrder[ 1 ],
+                                                              newOrder[ 2 ], newOrder[ 3 ],
+                                                              newOrder[ 4 ], newOrder[ 5 ],
+                                                              newOrder[ 6 ], newOrder[ 7 ],
+                                                              newOrder[ 8 ]));
               }
             }
             else  //////// polygon
@@ -5321,8 +4999,8 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
               }
             }
 
-            while ( srcElements.Length() < myLastCreatedElems.Length() )
-              srcElements.Append( *srcEdge );
+            while ( srcElements.size() < myLastCreatedElems.size() )
+              srcElements.push_back( *srcEdge );
 
           }  // loop on free faces
 
@@ -5359,8 +5037,8 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
 
         AddElement( nodeVec, anyFace.Init( elem ));
 
-        while ( srcElements.Length() < myLastCreatedElems.Length() )
-          srcElements.Append( elem );
+        while ( srcElements.size() < myLastCreatedElems.size() )
+          srcElements.push_back( elem );
       }
     }
   } // loop on swept elements
@@ -5380,13 +5058,17 @@ SMESH_MeshEditor::RotationSweep(TIDSortedElemSet   theElemSets[2],
                                 const bool         theMakeGroups,
                                 const bool         theMakeWalls)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
+
+  setElemsFirst( theElemSets );
+  myLastCreatedElems.reserve( theElemSets[0].size() * theNbSteps );
+  myLastCreatedNodes.reserve( theElemSets[1].size() * theNbSteps );
 
   // source elements for each generated one
   SMESH_SequenceOfElemPtr srcElems, srcNodes;
+  srcElems.reserve( theElemSets[0].size() );
+  srcNodes.reserve( theElemSets[1].size() );
 
-  MESSAGE( "RotationSweep()");
   gp_Trsf aTrsf;
   aTrsf.SetRotation( theAxis, theAngle );
   gp_Trsf aTrsf2;
@@ -5405,7 +5087,6 @@ SMESH_MeshEditor::RotationSweep(TIDSortedElemSet   theElemSets[2],
                                      myMesh->NbFaces(ORDER_QUADRATIC) +
                                      myMesh->NbVolumes(ORDER_QUADRATIC) );
   // loop on theElemSets
-  setElemsFirst( theElemSets );
   TIDSortedElemSet::iterator itElem;
   for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
   {
@@ -5457,8 +5138,8 @@ SMESH_MeshEditor::RotationSweep(TIDSortedElemSet   theElemSets[2],
               {
                 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
                 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
-                myLastCreatedNodes.Append(newNode);
-                srcNodes.Append( node );
+                myLastCreatedNodes.push_back(newNode);
+                srcNodes.push_back( node );
                 listNewNodes.push_back( newNode );
                 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
               }
@@ -5467,8 +5148,8 @@ SMESH_MeshEditor::RotationSweep(TIDSortedElemSet   theElemSets[2],
               }
               // create a corner node
               newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
-              myLastCreatedNodes.Append(newNode);
-              srcNodes.Append( node );
+              myLastCreatedNodes.push_back(newNode);
+              srcNodes.push_back( node );
               listNewNodes.push_back( newNode );
             }
             else {
@@ -5500,11 +5181,14 @@ SMESH_MeshEditor::RotationSweep(TIDSortedElemSet   theElemSets[2],
 //purpose  : standard construction
 //=======================================================================
 
-SMESH_MeshEditor::ExtrusParam::ExtrusParam( const gp_Vec&  theStep,
-                                            const int      theNbSteps,
-                                            const int      theFlags,
-                                            const double   theTolerance):
+SMESH_MeshEditor::ExtrusParam::ExtrusParam( const gp_Vec&            theStep,
+                                            const int                theNbSteps,
+                                            const std::list<double>& theScales,
+                                            const gp_XYZ*            theBasePoint,
+                                            const int                theFlags,
+                                            const double             theTolerance):
   myDir( theStep ),
+  myBaseP( Precision::Infinite(), 0, 0 ),
   myFlags( theFlags ),
   myTolerance( theTolerance ),
   myElemsToUse( NULL )
@@ -5514,6 +5198,37 @@ SMESH_MeshEditor::ExtrusParam::ExtrusParam( const gp_Vec&  theStep,
   for (int i=1; i<=theNbSteps; i++ )
     mySteps->Append( stepSize );
 
+  int nbScales = theScales.size();
+  if ( nbScales > 0 )
+  {
+    if ( IsLinearVariation() && nbScales < theNbSteps )
+    {
+      myScales.reserve( theNbSteps );
+      std::list<double>::const_iterator scale = theScales.begin();
+      double prevScale = 1.0;
+      for ( int iSc = 1; scale != theScales.end(); ++scale, ++iSc )
+      {
+        int      iStep = int( iSc / double( nbScales ) * theNbSteps + 0.5 );
+        int    stDelta = Max( 1, iStep - myScales.size());
+        double scDelta = ( *scale - prevScale ) / stDelta;
+        for ( int iStep = 0; iStep < stDelta; ++iStep )
+        {
+          myScales.push_back( prevScale + scDelta );
+          prevScale = myScales.back();
+        }
+        prevScale = *scale;
+      }
+    }
+    else
+    {
+      myScales.assign( theScales.begin(), theScales.end() );
+    }
+  }
+  if ( theBasePoint )
+  {
+    myBaseP = *theBasePoint;
+  }
+
   if (( theFlags & EXTRUSION_FLAG_SEW ) &&
       ( theTolerance > 0 ))
   {
@@ -5582,12 +5297,38 @@ SMESH_MeshEditor::ExtrusParam::ExtrusParam( const double theStepSize,
 //=======================================================================
 //function : ExtrusParam::SetElementsToUse
 //purpose  : stores elements to use for extrusion by normal, depending on
-//           state of EXTRUSION_FLAG_USE_INPUT_ELEMS_ONLY flag
+//           state of EXTRUSION_FLAG_USE_INPUT_ELEMS_ONLY flag;
+//           define myBaseP for scaling
 //=======================================================================
 
-void SMESH_MeshEditor::ExtrusParam::SetElementsToUse( const TIDSortedElemSet& elems )
+void SMESH_MeshEditor::ExtrusParam::SetElementsToUse( const TIDSortedElemSet& elems,
+                                                      const TIDSortedElemSet& nodes )
 {
   myElemsToUse = ToUseInpElemsOnly() ? & elems : 0;
+
+  if ( Precision::IsInfinite( myBaseP.X() )) // myBaseP not defined
+  {
+    myBaseP.SetCoord( 0.,0.,0. );
+    TIDSortedElemSet newNodes;
+
+    const TIDSortedElemSet* elemSets[] = { &elems, &nodes };
+    for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
+    {
+      const TIDSortedElemSet& elements = *( elemSets[ is2ndSet ]);
+      TIDSortedElemSet::const_iterator itElem = elements.begin();
+      for ( ; itElem != elements.end(); itElem++ )
+      {
+        const SMDS_MeshElement* elem = *itElem;
+        SMDS_ElemIteratorPtr     itN = elem->nodesIterator();
+        while ( itN->more() ) {
+          const SMDS_MeshElement* node = itN->next();
+          if ( newNodes.insert( node ).second )
+            myBaseP += SMESH_NodeXYZ( node );
+        }
+      }
+    }
+    myBaseP /= newNodes.size();
+  }
 }
 
 //=======================================================================
@@ -5648,7 +5389,7 @@ makeNodesByDir( SMESHDS_Mesh*                     mesh,
                 std::list<const SMDS_MeshNode*> & newNodes,
                 const bool                        makeMediumNodes)
 {
-  gp_XYZ p = SMESH_TNodeXYZ( srcNode );
+  gp_XYZ p = SMESH_NodeXYZ( srcNode );
 
   int nbNodes = 0;
   for ( beginStepIter( makeMediumNodes ); moreSteps(); ++nbNodes ) // loop on steps
@@ -5657,6 +5398,41 @@ makeNodesByDir( SMESHDS_Mesh*                     mesh,
     const SMDS_MeshNode * newNode = mesh->AddNode( p.X(), p.Y(), p.Z() );
     newNodes.push_back( newNode );
   }
+
+  if ( !myScales.empty() )
+  {
+    if ( makeMediumNodes && myMediumScales.empty() )
+    {
+      myMediumScales.resize( myScales.size() );
+      double prevFactor = 1.;
+      for ( size_t i = 0; i < myScales.size(); ++i )
+      {
+        myMediumScales[i] = 0.5 * ( prevFactor + myScales[i] );
+        prevFactor = myScales[i];
+      }
+    }
+    typedef std::vector<double>::iterator ScaleIt;
+    ScaleIt scales[] = { myScales.begin(), myMediumScales.begin() };
+
+    size_t iSc = 0, nbScales = myScales.size() + myMediumScales.size();
+
+    gp_XYZ center = myBaseP;
+    std::list<const SMDS_MeshNode*>::iterator nIt = newNodes.begin();
+    size_t iN  = 0;
+    for ( beginStepIter( makeMediumNodes ); moreSteps() && ( iN < nbScales ); ++nIt, ++iN )
+    {
+      center += myDir.XYZ() * nextStep();
+
+      iSc += int( makeMediumNodes );
+      ScaleIt& scale = scales[ iSc % 2 ];
+      
+      gp_XYZ xyz = SMESH_NodeXYZ( *nIt );
+      xyz = ( *scale * ( xyz - center )) + center;
+      mesh->MoveNode( *nIt, xyz.X(), xyz.Y(), xyz.Z() );
+
+      ++scale;
+    }
+  }
   return nbNodes;
 }
 
@@ -5671,7 +5447,7 @@ makeNodesByDirAndSew( SMESHDS_Mesh*                     mesh,
                       std::list<const SMDS_MeshNode*> & newNodes,
                       const bool                        makeMediumNodes)
 {
-  gp_XYZ P1 = SMESH_TNodeXYZ( srcNode );
+  gp_XYZ P1 = SMESH_NodeXYZ( srcNode );
 
   int nbNodes = 0;
   for ( beginStepIter( makeMediumNodes ); moreSteps(); ++nbNodes ) // loop on steps
@@ -5679,13 +5455,14 @@ makeNodesByDirAndSew( SMESHDS_Mesh*                     mesh,
     P1 += myDir.XYZ() * nextStep();
 
     // try to search in sequence of existing nodes
-    // if myNodes.Length()>0 we 'nave to use given sequence
+    // if myNodes.size()>0 we 'nave to use given sequence
     // else - use all nodes of mesh
     const SMDS_MeshNode * node = 0;
-    if ( myNodes.Length() > 0 ) {
-      int i;
-      for(i=1; i<=myNodes.Length(); i++) {
-        gp_XYZ P2 = SMESH_TNodeXYZ( myNodes.Value(i) );
+    if ( myNodes.Length() > 0 )
+    {
+      for ( int i = 1; i <= myNodes.Length(); i++ )
+      {
+        SMESH_NodeXYZ P2 = myNodes.Value(i);
         if (( P1 - P2 ).SquareModulus() < myTolerance * myTolerance )
         {
           node = myNodes.Value(i);
@@ -5693,10 +5470,12 @@ makeNodesByDirAndSew( SMESHDS_Mesh*                     mesh,
         }
       }
     }
-    else {
+    else
+    {
       SMDS_NodeIteratorPtr itn = mesh->nodesIterator();
-      while(itn->more()) {
-        SMESH_TNodeXYZ P2( itn->next() );
+      while(itn->more())
+      {
+        SMESH_NodeXYZ P2 = itn->next();
         if (( P1 - P2 ).SquareModulus() < myTolerance * myTolerance )
         {
           node = P2._node;
@@ -5728,7 +5507,7 @@ makeNodesByNormal2D( SMESHDS_Mesh*                     mesh,
 {
   const bool alongAvgNorm = ( myFlags & EXTRUSION_FLAG_BY_AVG_NORMAL );
 
-  gp_XYZ p = SMESH_TNodeXYZ( srcNode );
+  gp_XYZ p = SMESH_NodeXYZ( srcNode );
 
   // get normals to faces sharing srcNode
   vector< gp_XYZ > norms, baryCenters;
@@ -5748,7 +5527,7 @@ makeNodesByNormal2D( SMESHDS_Mesh*                     mesh,
         gp_XYZ bc(0,0,0);
         int nbN = 0;
         for ( SMDS_ElemIteratorPtr nIt = face->nodesIterator(); nIt->more(); ++nbN )
-          bc += SMESH_TNodeXYZ( nIt->next() );
+          bc += SMESH_NodeXYZ( nIt->next() );
         baryCenters.push_back( bc / nbN );
       }
     }
@@ -5831,7 +5610,7 @@ SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet     theElems[2],
                                   const int            theFlags,
                                   const double         theTolerance)
 {
-  ExtrusParam aParams( theStep, theNbSteps, theFlags, theTolerance );
+  ExtrusParam aParams( theStep, theNbSteps, std::list<double>(), 0, theFlags, theTolerance );
   return ExtrusionSweep( theElems, aParams, newElemsMap );
 }
 
@@ -5846,22 +5625,22 @@ SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet     theElemSets[2],
                                   ExtrusParam&         theParams,
                                   TTElemOfElemListMap& newElemsMap)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
+
+  setElemsFirst( theElemSets );
+  myLastCreatedElems.reserve( theElemSets[0].size() * theParams.NbSteps() );
+  myLastCreatedNodes.reserve( theElemSets[1].size() * theParams.NbSteps() );
 
   // source elements for each generated one
   SMESH_SequenceOfElemPtr srcElems, srcNodes;
+  srcElems.reserve( theElemSets[0].size() );
+  srcNodes.reserve( theElemSets[1].size() );
 
-  //SMESHDS_Mesh* aMesh = GetMeshDS();
-
-  setElemsFirst( theElemSets );
   const int nbSteps = theParams.NbSteps();
-  theParams.SetElementsToUse( theElemSets[0] );
+  theParams.SetElementsToUse( theElemSets[0], theElemSets[1] );
 
-  TNodeOfNodeListMap mapNewNodes;
-  //TNodeOfNodeVecMap mapNewNodes;
+  TNodeOfNodeListMap   mapNewNodes;
   TElemOfVecOfNnlmiMap mapElemNewNodes;
-  //TElemOfVecOfMapNodesMap mapElemNewNodes;
 
   const bool isQuadraticMesh = bool( myMesh->NbEdges(ORDER_QUADRATIC) +
                                      myMesh->NbFaces(ORDER_QUADRATIC) +
@@ -5915,8 +5694,8 @@ SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet     theElemSets[2],
             list<const SMDS_MeshNode*>::iterator newNodesIt = listNewNodes.begin();
             for ( ; newNodesIt != listNewNodes.end(); ++newNodesIt )
             {
-              myLastCreatedNodes.Append( *newNodesIt );
-              srcNodes.Append( node );
+              myLastCreatedNodes.push_back( *newNodesIt );
+              srcNodes.push_back( node );
             }
           }
           else
@@ -5957,9 +5736,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet     theElements[2],
                                        const gp_Pnt&        theRefPoint,
                                        const bool           theMakeGroups)
 {
-  MESSAGE("ExtrusionAlongTrack");
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
   int aNbE;
   std::list<double> aPrms;
@@ -6018,14 +5795,13 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet     theElements[2],
       return EXTR_BAD_STARTING_NODE;
     aItN = pSubMeshDS->GetNodes();
     while ( aItN->more() ) {
-      const SMDS_MeshNode* pNode = aItN->next();
-      const SMDS_EdgePosition* pEPos =
-        static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
+      const SMDS_MeshNode*  pNode = aItN->next();
+      SMDS_EdgePositionPtr pEPos = pNode->GetPosition();
       double aT = pEPos->GetUParameter();
       aPrms.push_back( aT );
     }
     //Extrusion_Error err =
-    MakeEdgePathPoints(aPrms, aTrackEdge, (aN1==theN1), fullList);
+    makeEdgePathPoints(aPrms, aTrackEdge, (aN1==theN1), fullList);
   } else if( aS.ShapeType() == TopAbs_WIRE ) {
     list< SMESH_subMesh* > LSM;
     TopTools_SequenceOfShape Edges;
@@ -6063,17 +5839,16 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet     theElements[2],
         aItN = locMeshDS->GetNodes();
         while ( aItN->more() ) {
           const SMDS_MeshNode* pNode = aItN->next();
-          const SMDS_EdgePosition* pEPos =
-            static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
+          SMDS_EdgePositionPtr pEPos = pNode->GetPosition();
           double aT = pEPos->GetUParameter();
           aPrms.push_back( aT );
         }
         list<SMESH_MeshEditor_PathPoint> LPP;
         //Extrusion_Error err =
-        MakeEdgePathPoints(aPrms, aTrackEdge,(aN1->GetID()==startNid), LPP);
+        makeEdgePathPoints(aPrms, aTrackEdge,(aN1->GetID()==startNid), LPP);
         LLPPs.push_back(LPP);
         UsedNums.Add(k);
-        // update startN for search following egde
+        // update startN for search following edge
         if( aN1->GetID() == startNid ) startNid = aN2->GetID();
         else startNid = aN1->GetID();
         break;
@@ -6113,7 +5888,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet     theElements[2],
     return EXTR_BAD_PATH_SHAPE;
   }
 
-  return MakeExtrElements(theElements, fullList, theHasAngles, theAngles, theLinearVariation,
+  return makeExtrElements(theElements, fullList, theHasAngles, theAngles, theLinearVariation,
                           theHasRefPoint, theRefPoint, theMakeGroups);
 }
 
@@ -6133,8 +5908,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet     theElements[2],
                                        const gp_Pnt&        theRefPoint,
                                        const bool           theMakeGroups)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
   int aNbE;
   std::list<double> aPrms;
@@ -6187,7 +5961,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet     theElements[2],
     SMDS_ElemIteratorPtr nIt;
 
     //check start node
-    if( !theTrack->GetMeshDS()->Contains(theN1) ) {
+    if( !theTrack->GetMeshDS()->Contains( theN1 )) {
       return EXTR_BAD_STARTING_NODE;
     }
 
@@ -6247,12 +6021,12 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet     theElements[2],
     int startNid = theN1->GetID();
     for ( size_t i = 1; i < aNodesList.size(); i++ )
     {
-      gp_Pnt     p1 = SMESH_TNodeXYZ( aNodesList[i-1] );
-      gp_Pnt     p2 = SMESH_TNodeXYZ( aNodesList[i] );
+      gp_Pnt     p1 = SMESH_NodeXYZ( aNodesList[i-1] );
+      gp_Pnt     p2 = SMESH_NodeXYZ( aNodesList[i] );
       TopoDS_Edge e = BRepBuilderAPI_MakeEdge( p1, p2 );
       list<SMESH_MeshEditor_PathPoint> LPP;
       aPrms.clear();
-      MakeEdgePathPoints(aPrms, e, (aNodesList[i-1]->GetID()==startNid), LPP);
+      makeEdgePathPoints(aPrms, e, (aNodesList[i-1]->GetID()==startNid), LPP);
       LLPPs.push_back(LPP);
       if ( aNodesList[i-1]->GetID() == startNid ) startNid = aNodesList[i  ]->GetID();
       else                                        startNid = aNodesList[i-1]->GetID();
@@ -6304,13 +6078,12 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet     theElements[2],
     while ( aItN->more() ) {
       const SMDS_MeshNode* pNode = aItN->next();
       if( pNode==aN1 || pNode==aN2 ) continue;
-      const SMDS_EdgePosition* pEPos =
-        static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
+      SMDS_EdgePositionPtr pEPos = pNode->GetPosition();
       double aT = pEPos->GetUParameter();
       aPrms.push_back( aT );
     }
     //Extrusion_Error err =
-    MakeEdgePathPoints(aPrms, aTrackEdge, (aN1==theN1), fullList);
+    makeEdgePathPoints(aPrms, aTrackEdge, (aN1==theN1), fullList);
   }
   else if( aS.ShapeType() == TopAbs_WIRE ) {
     list< SMESH_subMesh* > LSM;
@@ -6361,18 +6134,17 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet     theElements[2],
         aPrms.clear();
         aItN = locMeshDS->GetNodes();
         while ( aItN->more() ) {
-          const SMDS_MeshNode*     pNode = aItN->next();
-          const SMDS_EdgePosition* pEPos =
-            static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
+          const SMDS_MeshNode*  pNode = aItN->next();
+          SMDS_EdgePositionPtr pEPos = pNode->GetPosition();
           double aT = pEPos->GetUParameter();
           aPrms.push_back( aT );
         }
         list<SMESH_MeshEditor_PathPoint> LPP;
         //Extrusion_Error err =
-        MakeEdgePathPoints(aPrms, aTrackEdge, aN1isOK, LPP);
+        makeEdgePathPoints(aPrms, aTrackEdge, aN1isOK, LPP);
         LLPPs.push_back(LPP);
         UsedNums.Add(k);
-        // update startN for search following egde
+        // update startN for search following edge
         if ( aN1isOK ) aVprev = aV2;
         else           aVprev = aV1;
         break;
@@ -6405,17 +6177,17 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet     theElements[2],
     return EXTR_BAD_PATH_SHAPE;
   }
 
-  return MakeExtrElements(theElements, fullList, theHasAngles, theAngles, theLinearVariation,
+  return makeExtrElements(theElements, fullList, theHasAngles, theAngles, theLinearVariation,
                           theHasRefPoint, theRefPoint, theMakeGroups);
 }
 
 
 //=======================================================================
-//function : MakeEdgePathPoints
-//purpose  : auxilary for ExtrusionAlongTrack
+//function : makeEdgePathPoints
+//purpose  : auxiliary for ExtrusionAlongTrack
 //=======================================================================
 SMESH_MeshEditor::Extrusion_Error
-SMESH_MeshEditor::MakeEdgePathPoints(std::list<double>&                aPrms,
+SMESH_MeshEditor::makeEdgePathPoints(std::list<double>&                aPrms,
                                      const TopoDS_Edge&                aTrackEdge,
                                      bool                              FirstIsStart,
                                      list<SMESH_MeshEditor_PathPoint>& LPP)
@@ -6466,11 +6238,11 @@ SMESH_MeshEditor::MakeEdgePathPoints(std::list<double>&                aPrms,
 
 
 //=======================================================================
-//function : MakeExtrElements
-//purpose  : auxilary for ExtrusionAlongTrack
+//function : makeExtrElements
+//purpose  : auxiliary for ExtrusionAlongTrack
 //=======================================================================
 SMESH_MeshEditor::Extrusion_Error
-SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet                  theElemSets[2],
+SMESH_MeshEditor::makeExtrElements(TIDSortedElemSet                  theElemSets[2],
                                    list<SMESH_MeshEditor_PathPoint>& fullList,
                                    const bool                        theHasAngles,
                                    list<double>&                     theAngles,
@@ -6483,7 +6255,7 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet                  theElemSets
 
   // Angles
   if( theHasAngles && !theAngles.empty() && theLinearVariation )
-    LinearAngleVariation(aNbTP-1, theAngles);
+    linearAngleVariation(aNbTP-1, theAngles);
 
   // fill vector of path points with angles
   vector<SMESH_MeshEditor_PathPoint> aPPs;
@@ -6521,7 +6293,7 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet                  theElemSets
         while ( itN->more() ) {
           const SMDS_MeshElement* node = itN->next();
           if ( newNodes.insert( node ).second )
-            aGC += SMESH_TNodeXYZ( node );
+            aGC += SMESH_NodeXYZ( node );
         }
       }
     }
@@ -6565,7 +6337,7 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet                  theElemSets
           aTolAng=1.e-4;
 
           aV0x = aV0;
-          aPN0 = SMESH_TNodeXYZ( node );
+          aPN0 = SMESH_NodeXYZ( node );
 
           const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0];
           aP0x = aPP0.Pnt();
@@ -6582,7 +6354,7 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet                  theElemSets
             gp_Vec aV01x( aP0x, aP1x );
             aTrsf.SetTranslation( aV01x );
 
-            // traslated point
+            // translated point
             aV1x = aV0x.Transformed( aTrsf );
             aPN1 = aPN0.Transformed( aTrsf );
 
@@ -6613,13 +6385,13 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet                  theElemSets
               // create additional node
               gp_XYZ midP = 0.5 * ( aPN1.XYZ() + aPN0.XYZ() );
               const SMDS_MeshNode* newNode = aMesh->AddNode( midP.X(), midP.Y(), midP.Z() );
-              myLastCreatedNodes.Append(newNode);
-              srcNodes.Append( node );
+              myLastCreatedNodes.push_back(newNode);
+              srcNodes.push_back( node );
               listNewNodes.push_back( newNode );
             }
             const SMDS_MeshNode* newNode = aMesh->AddNode( aPN1.X(), aPN1.Y(), aPN1.Z() );
-            myLastCreatedNodes.Append(newNode);
-            srcNodes.Append( node );
+            myLastCreatedNodes.push_back(newNode);
+            srcNodes.push_back( node );
             listNewNodes.push_back( newNode );
 
             aPN0 = aPN1;
@@ -6645,8 +6417,8 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet                  theElemSets
               double y = ( N->Y() + P.Y() )/2.;
               double z = ( N->Z() + P.Z() )/2.;
               const SMDS_MeshNode* newN = aMesh->AddNode(x,y,z);
-              srcNodes.Append( node );
-              myLastCreatedNodes.Append(newN);
+              srcNodes.push_back( node );
+              myLastCreatedNodes.push_back(newN);
               aNodes[2*i] = newN;
               aNodes[2*i+1] = N;
               P = gp_XYZ(N->X(),N->Y(),N->Z());
@@ -6676,25 +6448,24 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet                  theElemSets
 
 
 //=======================================================================
-//function : LinearAngleVariation
-//purpose  : auxilary for ExtrusionAlongTrack
+//function : linearAngleVariation
+//purpose  : spread values over nbSteps
 //=======================================================================
-void SMESH_MeshEditor::LinearAngleVariation(const int nbSteps,
+
+void SMESH_MeshEditor::linearAngleVariation(const int     nbSteps,
                                             list<double>& Angles)
 {
   int nbAngles = Angles.size();
-  if( nbSteps > nbAngles ) {
+  if( nbSteps > nbAngles && nbAngles > 0 )
+  {
     vector<double> theAngles(nbAngles);
-    list<double>::iterator it = Angles.begin();
-    int i = -1;
-    for(; it!=Angles.end(); it++) {
-      i++;
-      theAngles[i] = (*it);
-    }
+    theAngles.assign( Angles.begin(), Angles.end() );
+
     list<double> res;
     double rAn2St = double( nbAngles ) / double( nbSteps );
     double angPrev = 0, angle;
-    for ( int iSt = 0; iSt < nbSteps; ++iSt ) {
+    for ( int iSt = 0; iSt < nbSteps; ++iSt )
+    {
       double angCur = rAn2St * ( iSt+1 );
       double angCurFloor  = floor( angCur );
       double angPrevFloor = floor( angPrev );
@@ -6716,10 +6487,7 @@ void SMESH_MeshEditor::LinearAngleVariation(const int nbSteps,
       res.push_back(angle);
       angPrev = angCur;
     }
-    Angles.clear();
-    it = res.begin();
-    for(; it!=res.end(); it++)
-      Angles.push_back( *it );
+    Angles.swap( res );
   }
 }
 
@@ -6743,44 +6511,36 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
                              const bool         theMakeGroups,
                              SMESH_Mesh*        theTargetMesh)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
+  myLastCreatedElems.reserve( theElems.size() );
 
   bool needReverse = false;
   string groupPostfix;
   switch ( theTrsf.Form() ) {
   case gp_PntMirror:
-    MESSAGE("gp_PntMirror");
     needReverse = true;
     groupPostfix = "mirrored";
     break;
   case gp_Ax1Mirror:
-    MESSAGE("gp_Ax1Mirror");
     groupPostfix = "mirrored";
     break;
   case gp_Ax2Mirror:
-    MESSAGE("gp_Ax2Mirror");
     needReverse = true;
     groupPostfix = "mirrored";
     break;
   case gp_Rotation:
-    MESSAGE("gp_Rotation");
     groupPostfix = "rotated";
     break;
   case gp_Translation:
-    MESSAGE("gp_Translation");
     groupPostfix = "translated";
     break;
   case gp_Scale:
-    MESSAGE("gp_Scale");
     groupPostfix = "scaled";
     break;
   case gp_CompoundTrsf: // different scale by axis
-    MESSAGE("gp_CompoundTrsf");
     groupPostfix = "scaled";
     break;
   default:
-    MESSAGE("default");
     needReverse = false;
     groupPostfix = "transformed";
   }
@@ -6847,14 +6607,14 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
         if ( theTargetMesh ) {
           const SMDS_MeshNode * newNode = aTgtMesh->AddNode( coord[0], coord[1], coord[2] );
           n2n_isnew.first->second = newNode;
-          myLastCreatedNodes.Append(newNode);
-          srcNodes.Append( node );
+          myLastCreatedNodes.push_back(newNode);
+          srcNodes.push_back( node );
         }
         else if ( theCopy ) {
           const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
           n2n_isnew.first->second = newNode;
-          myLastCreatedNodes.Append(newNode);
-          srcNodes.Append( node );
+          myLastCreatedNodes.push_back(newNode);
+          srcNodes.push_back( node );
         }
         else {
           aMesh->MoveNode( node, coord[0], coord[1], coord[2] );
@@ -6897,8 +6657,8 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
 
     if ( geomType == SMDSGeom_POLYHEDRA )  // ------------------ polyhedral volume
     {
-      const SMDS_VtkVolume* aPolyedre = dynamic_cast<const SMDS_VtkVolume*>( elem );
-      if (!aPolyedre)
+      const SMDS_MeshVolume* aPolyedre = SMDS_Mesh::DownCast< SMDS_MeshVolume >( elem );
+      if ( !aPolyedre )
         continue;
       nodes.clear();
       bool allTransformed = true;
@@ -6944,7 +6704,7 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
     if ( editor ) {
       // copy in this or a new mesh
       if ( editor->AddElement( nodes, elemType.Init( elem, /*basicOnly=*/false )))
-        srcElems.Append( elem );
+        srcElems.push_back( elem );
     }
     else {
       // reverse element as it was reversed by transformation
@@ -6955,7 +6715,7 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
   } // loop on elements
 
   if ( editor && editor != this )
-    myLastCreatedElems = editor->myLastCreatedElems;
+    myLastCreatedElems.swap( editor->myLastCreatedElems );
 
   PGroupIDs newGroupIDs;
 
@@ -6966,14 +6726,110 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
   return newGroupIDs;
 }
 
-//=======================================================================
+//================================================================================
+/*!
+ * \brief Make an offset mesh from a source 2D mesh
+ *  \param [in] theElements - source faces
+ *  \param [in] theValue - offset value
+ *  \param [out] theTgtMesh - a mesh to add offset elements to
+ *  \param [in] theMakeGroups - to generate groups
+ *  \return PGroupIDs - IDs of created groups
+ */
+//================================================================================
+
+SMESH_MeshEditor::PGroupIDs SMESH_MeshEditor::Offset( TIDSortedElemSet & theElements,
+                                                      const double       theValue,
+                                                      SMESH_Mesh*        theTgtMesh,
+                                                      const bool         theMakeGroups,
+                                                      const bool         theCopyElements,
+                                                      const bool         theFixSelfIntersection)
+{
+  SMESHDS_Mesh*    meshDS = GetMeshDS();
+  SMESHDS_Mesh* tgtMeshDS = theTgtMesh->GetMeshDS();
+  SMESH_MeshEditor tgtEditor( theTgtMesh );
+
+  SMDS_ElemIteratorPtr eIt;
+  if ( theElements.empty() ) eIt = meshDS->elementsIterator( SMDSAbs_Face );
+  else                       eIt = SMESHUtils::elemSetIterator( theElements );
+
+  SMESH_MeshAlgos::TEPairVec new2OldFaces;
+  SMESH_MeshAlgos::TNPairVec new2OldNodes;
+  std::unique_ptr< SMDS_Mesh > offsetMesh
+    ( SMESH_MeshAlgos::MakeOffset( eIt, *meshDS, theValue,
+                                   theFixSelfIntersection,
+                                   new2OldFaces, new2OldNodes ));
+  if ( offsetMesh->NbElements() == 0 )
+    return PGroupIDs(); // MakeOffset() failed
+
+
+  if ( theTgtMesh == myMesh && !theCopyElements )
+  {
+    // clear the source elements
+    if ( theElements.empty() ) eIt = meshDS->elementsIterator( SMDSAbs_Face );
+    else                       eIt = SMESHUtils::elemSetIterator( theElements );
+    while ( eIt->more() )
+      meshDS->RemoveFreeElement( eIt->next(), 0 );
+  }
+
+  // offsetMesh->Modified();
+  // offsetMesh->CompactMesh(); // make IDs start from 1
+
+  // source elements for each generated one
+  SMESH_SequenceOfElemPtr srcElems, srcNodes;
+  srcElems.reserve( new2OldFaces.size() );
+  srcNodes.reserve( new2OldNodes.size() );
+
+  ClearLastCreated();
+  myLastCreatedElems.reserve( new2OldFaces.size() );
+  myLastCreatedNodes.reserve( new2OldNodes.size() );
+
+  // copy offsetMesh to theTgtMesh
+
+  int idShift = meshDS->MaxNodeID();
+  for ( size_t i = 0; i < new2OldNodes.size(); ++i )
+    if ( const SMDS_MeshNode* n = new2OldNodes[ i ].first )
+    {
+      if ( n->NbInverseElements() > 0 )
+      {
+        const SMDS_MeshNode* n2 =
+          tgtMeshDS->AddNodeWithID( n->X(), n->Y(), n->Z(), idShift + n->GetID() );
+        myLastCreatedNodes.push_back( n2 );
+        srcNodes.push_back( new2OldNodes[ i ].second );
+      }
+    }
+
+  ElemFeatures elemType;
+  for ( size_t i = 0; i < new2OldFaces.size(); ++i )
+    if ( const SMDS_MeshElement* f = new2OldFaces[ i ].first )
+    {
+      elemType.Init( f );
+      elemType.myNodes.clear();
+      for ( SMDS_NodeIteratorPtr nIt = f->nodeIterator(); nIt->more(); )
+      {
+        const SMDS_MeshNode* n2 = nIt->next();
+        elemType.myNodes.push_back( tgtMeshDS->FindNode( idShift + n2->GetID() ));
+      }
+      tgtEditor.AddElement( elemType.myNodes, elemType );
+      srcElems.push_back( new2OldFaces[ i ].second );
+    }
+
+  myLastCreatedElems.swap( tgtEditor.myLastCreatedElems );
+
+  PGroupIDs newGroupIDs;
+  if ( theMakeGroups )
+    newGroupIDs = generateGroups( srcNodes, srcElems, "offset", theTgtMesh, false );
+
+  return newGroupIDs;
+}
+
+//=======================================================================
 /*!
  * \brief Create groups of elements made during transformation
  *  \param nodeGens - nodes making corresponding myLastCreatedNodes
  *  \param elemGens - elements making corresponding myLastCreatedElems
- *  \param postfix - to append to names of new groups
+ *  \param postfix - to push_back to names of new groups
  *  \param targetMesh - mesh to create groups in
- *  \param topPresent - is there "top" elements that are created by sweeping
+ *  \param topPresent - is there are "top" elements that are created by sweeping
  */
 //=======================================================================
 
@@ -7026,30 +6882,30 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
   {
     const SMESH_SequenceOfElemPtr& gens  = isNodes ? nodeGens : elemGens;
     const SMESH_SequenceOfElemPtr& elems = isNodes ? myLastCreatedNodes : myLastCreatedElems;
-    if ( gens.Length() != elems.Length() )
+    if ( gens.size() != elems.size() )
       throw SALOME_Exception("SMESH_MeshEditor::generateGroups(): invalid args");
 
     // loop on created elements
-    for (int iElem = 1; iElem <= elems.Length(); ++iElem )
+    for (size_t iElem = 0; iElem < elems.size(); ++iElem )
     {
-      const SMDS_MeshElement* sourceElem = gens( iElem );
+      const SMDS_MeshElement* sourceElem = gens[ iElem ];
       if ( !sourceElem ) {
         MESSAGE("generateGroups(): NULL source element");
         continue;
       }
       list< TOldNewGroup > & groupsOldNew = groupsByType[ sourceElem->GetType() ];
       if ( groupsOldNew.empty() ) { // no groups of this type at all
-        while ( iElem < gens.Length() && gens( iElem+1 ) == sourceElem )
+        while ( iElem+1 < gens.size() && gens[ iElem+1 ] == sourceElem )
           ++iElem; // skip all elements made by sourceElem
         continue;
       }
       // collect all elements made by the iElem-th sourceElem
       resultElems.clear();
-      if ( const SMDS_MeshElement* resElem = elems( iElem ))
+      if ( const SMDS_MeshElement* resElem = elems[ iElem ])
         if ( resElem != sourceElem )
           resultElems.push_back( resElem );
-      while ( iElem < gens.Length() && gens( iElem+1 ) == sourceElem )
-        if ( const SMDS_MeshElement* resElem = elems( ++iElem ))
+      while ( iElem+1 < gens.size() && gens[ iElem+1 ] == sourceElem )
+        if ( const SMDS_MeshElement* resElem = elems[ ++iElem ])
           if ( resElem != sourceElem )
             resultElems.push_back( resElem );
 
@@ -7089,7 +6945,7 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
           {
             SMDS_MeshGroup & newTopGroup = gOldNew->get<2>()->SMDSGroup();
             newTopGroup.Add( topElem );
-         }
+          }
         }
       }
     } // loop on created elements
@@ -7176,8 +7032,7 @@ void SMESH_MeshEditor::FindCoincidentNodes (TIDSortedNodeSet &   theNodes,
                                             TListOfListOfNodes & theGroupsOfNodes,
                                             bool                 theSeparateCornersAndMedium)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
   if ( myMesh->NbEdges  ( ORDER_QUADRATIC ) +
        myMesh->NbFaces  ( ORDER_QUADRATIC ) +
@@ -7190,7 +7045,7 @@ void SMESH_MeshEditor::FindCoincidentNodes (TIDSortedNodeSet &   theNodes,
   if ( theNodes.empty() ) // get all nodes in the mesh
   {
     TIDSortedNodeSet* nodes[2] = { &corners, &medium };
-    SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator(/*idInceasingOrder=*/true);
+    SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator();
     if ( theSeparateCornersAndMedium )
       while ( nIt->more() )
       {
@@ -7200,7 +7055,7 @@ void SMESH_MeshEditor::FindCoincidentNodes (TIDSortedNodeSet &   theNodes,
       }
     else
       while ( nIt->more() )
-        theNodes.insert( theNodes.end(),nIt->next() );
+        theNodes.insert( theNodes.end(), nIt->next() );
   }
   else if ( theSeparateCornersAndMedium ) // separate corners from medium nodes
   {
@@ -7233,76 +7088,52 @@ int SMESH_MeshEditor::SimplifyFace (const vector<const SMDS_MeshNode *>& faceNod
                                     vector<int>&                         quantities) const
 {
   int nbNodes = faceNodes.size();
-
-  if (nbNodes < 3)
+  while ( faceNodes[ 0 ] == faceNodes[ nbNodes-1 ] && nbNodes > 2 )
+    --nbNodes;
+  if ( nbNodes < 3 )
     return 0;
+  size_t prevNbQuant = quantities.size();
 
-  set<const SMDS_MeshNode*> nodeSet;
-
-  // get simple seq of nodes
-  vector<const SMDS_MeshNode*> simpleNodes( nbNodes );
-  int iSimple = 0;
+  vector< const SMDS_MeshNode* > simpleNodes; simpleNodes.reserve( nbNodes );
+  map< const SMDS_MeshNode*, int > nodeIndices; // indices within simpleNodes
+  map< const SMDS_MeshNode*, int >::iterator nInd;
 
-  simpleNodes[iSimple++] = faceNodes[0];
-  for (int iCur = 1; iCur < nbNodes; iCur++) {
-    if (faceNodes[iCur] != simpleNodes[iSimple - 1]) {
-      simpleNodes[iSimple++] = faceNodes[iCur];
-      nodeSet.insert( faceNodes[iCur] );
-    }
-  }
-  int nbUnique = nodeSet.size();
-  int nbSimple = iSimple;
-  if (simpleNodes[nbSimple - 1] == simpleNodes[0]) {
-    nbSimple--;
-    iSimple--;
-  }
-
-  if (nbUnique < 3)
-    return 0;
-
-  // separate loops
-  int nbNew = 0;
-  bool foundLoop = (nbSimple > nbUnique);
-  while (foundLoop) {
-    foundLoop = false;
-    set<const SMDS_MeshNode*> loopSet;
-    for (iSimple = 0; iSimple < nbSimple && !foundLoop; iSimple++) {
-      const SMDS_MeshNode* n = simpleNodes[iSimple];
-      if (!loopSet.insert( n ).second) {
-        foundLoop = true;
-
-        // separate loop
-        int iC = 0, curLast = iSimple;
-        for (; iC < curLast; iC++) {
-          if (simpleNodes[iC] == n) break;
-        }
-        int loopLen = curLast - iC;
-        if (loopLen > 2) {
-          // create sub-element
-          nbNew++;
-          quantities.push_back(loopLen);
-          for (; iC < curLast; iC++) {
-            poly_nodes.push_back(simpleNodes[iC]);
-          }
-        }
-        // shift the rest nodes (place from the first loop position)
-        for (iC = curLast + 1; iC < nbSimple; iC++) {
-          simpleNodes[iC - loopLen] = simpleNodes[iC];
+  nodeIndices.insert( make_pair( faceNodes[0], 0 ));
+  simpleNodes.push_back( faceNodes[0] );
+  for ( int iCur = 1; iCur < nbNodes; iCur++ )
+  {
+    if ( faceNodes[ iCur ] != simpleNodes.back() )
+    {
+      int index = simpleNodes.size();
+      nInd = nodeIndices.insert( make_pair( faceNodes[ iCur ], index )).first;
+      int prevIndex = nInd->second;
+      if ( prevIndex < index )
+      {
+        // a sub-loop found
+        int loopLen = index - prevIndex;
+        if ( loopLen > 2 )
+        {
+          // store the sub-loop
+          quantities.push_back( loopLen );
+          for ( int i = prevIndex; i < index; i++ )
+            poly_nodes.push_back( simpleNodes[ i ]);
         }
-        nbSimple -= loopLen;
-        iSimple -= loopLen;
+        simpleNodes.resize( prevIndex+1 );
       }
-    } // for (iSimple = 0; iSimple < nbSimple; iSimple++)
-  } // while (foundLoop)
+      else
+      {
+        simpleNodes.push_back( faceNodes[ iCur ]);
+      }
+    }
+  }
 
-  if (iSimple > 2) {
-    nbNew++;
-    quantities.push_back(iSimple);
-    for (int i = 0; i < iSimple; i++)
-      poly_nodes.push_back(simpleNodes[i]);
+  if ( simpleNodes.size() > 2 )
+  {
+    quantities.push_back( simpleNodes.size() );
+    poly_nodes.insert ( poly_nodes.end(), simpleNodes.begin(), simpleNodes.end() );
   }
 
-  return nbNew;
+  return quantities.size() - prevNbQuant;
 }
 
 //=======================================================================
@@ -7311,17 +7142,17 @@ int SMESH_MeshEditor::SimplifyFace (const vector<const SMDS_MeshNode *>& faceNod
 //           in all elements.
 //=======================================================================
 
-void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
+void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes,
+                                   const bool           theAvoidMakingHoles)
 {
-  MESSAGE("MergeNodes");
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
-  SMESHDS_Mesh* aMesh = GetMeshDS();
+  SMESHDS_Mesh* mesh = GetMeshDS();
 
   TNodeNodeMap nodeNodeMap; // node to replace - new node
   set<const SMDS_MeshElement*> elems; // all elements with changed nodes
   list< int > rmElemIds, rmNodeIds;
+  vector< ElemFeatures > newElemDefs;
 
   // Fill nodeNodeMap and elems
 
@@ -7335,17 +7166,6 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
     {
       const SMDS_MeshNode* nToRemove = *nIt;
       nodeNodeMap.insert( make_pair( nToRemove, nToKeep ));
-      if ( nToRemove != nToKeep )
-      {
-        rmNodeIds.push_back( nToRemove->GetID() );
-        AddToSameGroups( nToKeep, nToRemove, aMesh );
-        // set _alwaysComputed to a sub-mesh of VERTEX to enable mesh computing
-        // after MergeNodes() w/o creating node in place of merged ones.
-        const SMDS_PositionPtr& pos = nToRemove->GetPosition();
-        if ( pos && pos->GetTypeOfPosition() == SMDS_TOP_VERTEX )
-          if ( SMESH_subMesh* sm = myMesh->GetSubMeshContaining( nToRemove->getshapeId() ))
-            sm->SetIsAlwaysComputed( true );
-      }
       SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
       while ( invElemIt->more() ) {
         const SMDS_MeshElement* elem = invElemIt->next();
@@ -7353,611 +7173,624 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
       }
     }
   }
-  // Change element nodes or remove an element
 
-  set<const SMDS_MeshNode*> nodeSet;
-  vector< const SMDS_MeshNode*> curNodes, uniqueNodes;
-  vector<int> iRepl;
-  ElemFeatures elemType;
+  // Apply recursive replacements (BUG 0020185)
+  TNodeNodeMap::iterator nnIt = nodeNodeMap.begin();
+  for ( ; nnIt != nodeNodeMap.end(); ++nnIt )
+  {
+    const SMDS_MeshNode* nToKeep = nnIt->second;
+    TNodeNodeMap::iterator nnIt_i = nodeNodeMap.find( nToKeep );
+    while ( nnIt_i != nodeNodeMap.end() && nnIt_i->second != nnIt->second )
+    {
+      nToKeep = nnIt_i->second;
+      nnIt->second = nToKeep;
+      nnIt_i = nodeNodeMap.find( nToKeep );
+    }
+  }
+
+  if ( theAvoidMakingHoles )
+  {
+    // find elements whose topology changes
+
+    vector<const SMDS_MeshElement*> pbElems;
+    set<const SMDS_MeshElement*>::iterator eIt = elems.begin();
+    for ( ; eIt != elems.end(); ++eIt )
+    {
+      const SMDS_MeshElement* elem = *eIt;
+      SMDS_ElemIteratorPtr     itN = elem->nodesIterator();
+      while ( itN->more() )
+      {
+        const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( itN->next() );
+        TNodeNodeMap::iterator nnIt = nodeNodeMap.find( n );
+        if ( nnIt != nodeNodeMap.end() && elem->GetNodeIndex( nnIt->second ) >= 0 )
+        {
+          // several nodes of elem stick
+          pbElems.push_back( elem );
+          break;
+        }
+      }
+    }
+    // exclude from merge nodes causing spoiling element
+    for ( size_t iLoop = 0; iLoop < pbElems.size(); ++iLoop ) // avoid infinite cycle
+    {
+      bool nodesExcluded = false;
+      for ( size_t i = 0; i < pbElems.size(); ++i )
+      {
+        size_t prevNbMergeNodes = nodeNodeMap.size();
+        if ( !applyMerge( pbElems[i], newElemDefs, nodeNodeMap, /*noHoles=*/true ) &&
+             prevNbMergeNodes < nodeNodeMap.size() )
+          nodesExcluded = true;
+      }
+      if ( !nodesExcluded )
+        break;
+    }
+  }
+
+  for ( nnIt = nodeNodeMap.begin(); nnIt != nodeNodeMap.end(); ++nnIt )
+  {
+    const SMDS_MeshNode* nToRemove = nnIt->first;
+    const SMDS_MeshNode* nToKeep   = nnIt->second;
+    if ( nToRemove != nToKeep )
+    {
+      rmNodeIds.push_back( nToRemove->GetID() );
+      AddToSameGroups( nToKeep, nToRemove, mesh );
+      // set _alwaysComputed to a sub-mesh of VERTEX to enable further mesh computing
+      // w/o creating node in place of merged ones.
+      SMDS_PositionPtr pos = nToRemove->GetPosition();
+      if ( pos && pos->GetTypeOfPosition() == SMDS_TOP_VERTEX )
+        if ( SMESH_subMesh* sm = myMesh->GetSubMeshContaining( nToRemove->getshapeId() ))
+          sm->SetIsAlwaysComputed( true );
+    }
+  }
+
+  // Change element nodes or remove an element
 
   set<const SMDS_MeshElement*>::iterator eIt = elems.begin();
   for ( ; eIt != elems.end(); eIt++ )
   {
     const SMDS_MeshElement* elem = *eIt;
-    const           int  nbNodes = elem->NbNodes();
-    const           int aShapeId = FindShape( elem );
+    SMESHDS_SubMesh*          sm = mesh->MeshElements( elem->getshapeId() );
 
-    nodeSet.clear();
-    curNodes.resize( nbNodes );
-    uniqueNodes.resize( nbNodes );
-    iRepl.resize( nbNodes );
-    int iUnique = 0, iCur = 0, nbRepl = 0;
+    bool keepElem = applyMerge( elem, newElemDefs, nodeNodeMap, /*noHoles=*/false );
+    if ( !keepElem )
+      rmElemIds.push_back( elem->GetID() );
 
-    // get new seq of nodes
-    SMDS_ElemIteratorPtr itN = elem->nodesIterator();
-    while ( itN->more() )
-    {
-      const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( itN->next() );
-
-      TNodeNodeMap::iterator nnIt = nodeNodeMap.find( n );
-      if ( nnIt != nodeNodeMap.end() ) { // n sticks
-        n = (*nnIt).second;
-        { ////////// BUG 0020185: begin
-          bool stopRecur = false;
-          set<const SMDS_MeshNode*> nodesRecur;
-          nodesRecur.insert(n);
-          while (!stopRecur) {
-            TNodeNodeMap::iterator nnIt_i = nodeNodeMap.find( n );
-            if ( nnIt_i != nodeNodeMap.end() ) { // n sticks
-              n = (*nnIt_i).second;
-              if (!nodesRecur.insert(n).second) {
-                // error: recursive dependancy
-                stopRecur = true;
-              }
-            }
-            else
-              stopRecur = true;
-          }
-        } ////////// BUG 0020185: end
+    for ( size_t i = 0; i < newElemDefs.size(); ++i )
+    {
+      if ( i > 0 || !mesh->ChangeElementNodes( elem,
+                                               & newElemDefs[i].myNodes[0],
+                                               newElemDefs[i].myNodes.size() ))
+      {
+        if ( i == 0 )
+        {
+          newElemDefs[i].SetID( elem->GetID() );
+          mesh->RemoveFreeElement(elem, sm, /*fromGroups=*/false);
+          if ( !keepElem ) rmElemIds.pop_back();
+        }
+        else
+        {
+          newElemDefs[i].SetID( -1 );
+        }
+        SMDS_MeshElement* newElem = this->AddElement( newElemDefs[i].myNodes, newElemDefs[i] );
+        if ( sm && newElem )
+          sm->AddElement( newElem );
+        if ( elem != newElem )
+          ReplaceElemInGroups( elem, newElem, mesh );
       }
-      curNodes[ iCur ] = n;
-      bool isUnique = nodeSet.insert( n ).second;
-      if ( isUnique )
-        uniqueNodes[ iUnique++ ] = n;
-      else
-        iRepl[ nbRepl++ ] = iCur;
-      iCur++;
     }
+  }
+
+  // Remove bad elements, then equal nodes (order important)
+  Remove( rmElemIds, /*isNodes=*/false );
+  Remove( rmNodeIds, /*isNodes=*/true );
+
+  return;
+}
+
+//=======================================================================
+//function : applyMerge
+//purpose  : Compute new connectivity of an element after merging nodes
+//  \param [in] elems - the element
+//  \param [out] newElemDefs - definition(s) of result element(s)
+//  \param [inout] nodeNodeMap - nodes to merge
+//  \param [in] avoidMakingHoles - if true and and the element becomes invalid
+//              after merging (but not degenerated), removes nodes causing
+//              the invalidity from \a nodeNodeMap.
+//  \return bool - true if the element should be removed
+//=======================================================================
+
+bool SMESH_MeshEditor::applyMerge( const SMDS_MeshElement* elem,
+                                   vector< ElemFeatures >& newElemDefs,
+                                   TNodeNodeMap&           nodeNodeMap,
+                                   const bool              avoidMakingHoles )
+{
+  bool toRemove = false; // to remove elem
+  int nbResElems = 1;    // nb new elements
+
+  newElemDefs.resize(nbResElems);
+  newElemDefs[0].Init( elem );
+  newElemDefs[0].myNodes.clear();
+
+  set<const SMDS_MeshNode*> nodeSet;
+  vector< const SMDS_MeshNode*>   curNodes;
+  vector< const SMDS_MeshNode*> & uniqueNodes = newElemDefs[0].myNodes;
+  vector<int> iRepl;
+
+  const        int  nbNodes = elem->NbNodes();
+  SMDSAbs_EntityType entity = elem->GetEntityType();
+
+  curNodes.resize( nbNodes );
+  uniqueNodes.resize( nbNodes );
+  iRepl.resize( nbNodes );
+  int iUnique = 0, iCur = 0, nbRepl = 0;
+
+  // Get new seq of nodes
+
+  SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+  while ( itN->more() )
+  {
+    const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( itN->next() );
+
+    TNodeNodeMap::iterator nnIt = nodeNodeMap.find( n );
+    if ( nnIt != nodeNodeMap.end() ) {
+      n = (*nnIt).second;
+    }
+    curNodes[ iCur ] = n;
+    bool isUnique = nodeSet.insert( n ).second;
+    if ( isUnique )
+      uniqueNodes[ iUnique++ ] = n;
+    else
+      iRepl[ nbRepl++ ] = iCur;
+    iCur++;
+  }
 
-    // Analyse element topology after replacement
+  // Analyse element topology after replacement
 
-    bool isOk = true;
-    int nbUniqueNodes = nodeSet.size();
-    if ( nbNodes != nbUniqueNodes ) // some nodes stick
+  int nbUniqueNodes = nodeSet.size();
+  if ( nbNodes != nbUniqueNodes ) // some nodes stick
+  {
+    toRemove = true;
+    nbResElems = 0;
+
+    if ( newElemDefs[0].myIsQuad && newElemDefs[0].myType == SMDSAbs_Face && nbNodes > 6 )
     {
-      if (elem->IsPoly()) // Polygons and Polyhedral volumes
+      // if corner nodes stick, remove medium nodes between them from uniqueNodes
+      int nbCorners = nbNodes / 2;
+      for ( int iCur = 0; iCur < nbCorners; ++iCur )
       {
-        if (elem->GetType() == SMDSAbs_Face) // Polygon
+        int iNext = ( iCur + 1 ) % nbCorners;
+        if ( curNodes[ iCur ] == curNodes[ iNext ] ) // corners stick
         {
-          elemType.Init( elem );
-          const bool isQuad = elemType.myIsQuad;
-          if ( isQuad )
-            SMDS_MeshCell::applyInterlace // interlace medium and corner nodes
-              ( SMDS_MeshCell::interlacedSmdsOrder( SMDSEntity_Quad_Polygon, nbNodes ), curNodes );
-
-          // a polygon can divide into several elements
-          vector<const SMDS_MeshNode *> polygons_nodes;
-          vector<int> quantities;
-          int nbNew = SimplifyFace( curNodes, polygons_nodes, quantities );
-          if (nbNew > 0)
+          int iMedium = iCur + nbCorners;
+          vector< const SMDS_MeshNode* >::iterator i =
+            std::find( uniqueNodes.begin() + nbCorners - nbRepl,
+                       uniqueNodes.end(),
+                       curNodes[ iMedium ]);
+          if ( i != uniqueNodes.end() )
           {
-            vector<const SMDS_MeshNode *> face_nodes;
-            int inode = 0;
-            for (int iface = 0; iface < nbNew; iface++)
-            {
-              int nbNewNodes = quantities[iface];
-              face_nodes.assign( polygons_nodes.begin() + inode,
-                                 polygons_nodes.begin() + inode + nbNewNodes );
-              inode += nbNewNodes;
-              if ( isQuad ) // check if a result elem is a valid quadratic polygon
-              {
-                bool isValid = ( nbNewNodes % 2 == 0 );
-                for ( int i = 0; i < nbNewNodes && isValid; ++i )
-                  isValid = ( elem->IsMediumNode( face_nodes[i]) == bool( i % 2 ));
-                elemType.SetQuad( isValid );
-                if ( isValid ) // put medium nodes after corners
-                  SMDS_MeshCell::applyInterlaceRev
-                    ( SMDS_MeshCell::interlacedSmdsOrder( SMDSEntity_Quad_Polygon,
-                                                          nbNewNodes ), face_nodes );
-              }
-              elemType.SetPoly(( nbNewNodes / ( elemType.myIsQuad + 1 ) > 4 ));
-
-              SMDS_MeshElement* newElem = AddElement( face_nodes, elemType );
-              if ( aShapeId )
-                aMesh->SetMeshElementOnShape(newElem, aShapeId);
-            }
+            --nbUniqueNodes;
+            for ( ; i+1 != uniqueNodes.end(); ++i )
+              *i = *(i+1);
           }
-          rmElemIds.push_back(elem->GetID());
+        }
+      }
+    }
 
-        } // Polygon
+    switch ( entity )
+    {
+    case SMDSEntity_Polygon:
+    case SMDSEntity_Quad_Polygon: // Polygon
+    {
+      ElemFeatures* elemType = & newElemDefs[0];
+      const bool isQuad = elemType->myIsQuad;
+      if ( isQuad )
+        SMDS_MeshCell::applyInterlace // interlace medium and corner nodes
+          ( SMDS_MeshCell::interlacedSmdsOrder( SMDSEntity_Quad_Polygon, nbNodes ), curNodes );
 
-        else if (elem->GetType() == SMDSAbs_Volume) // Polyhedral volume
+      // a polygon can divide into several elements
+      vector<const SMDS_MeshNode *> polygons_nodes;
+      vector<int> quantities;
+      nbResElems = SimplifyFace( curNodes, polygons_nodes, quantities );
+      newElemDefs.resize( nbResElems );
+      for ( int inode = 0, iface = 0; iface < nbResElems; iface++ )
+      {
+        ElemFeatures* elemType = & newElemDefs[iface];
+        if ( iface ) elemType->Init( elem );
+
+        vector<const SMDS_MeshNode *>& face_nodes = elemType->myNodes;
+        int nbNewNodes = quantities[iface];
+        face_nodes.assign( polygons_nodes.begin() + inode,
+                           polygons_nodes.begin() + inode + nbNewNodes );
+        inode += nbNewNodes;
+        if ( isQuad ) // check if a result elem is a valid quadratic polygon
         {
-          if (nbUniqueNodes < 4) {
-            rmElemIds.push_back(elem->GetID());
-          }
-          else {
-            // each face has to be analyzed in order to check volume validity
-            const SMDS_VtkVolume* aPolyedre = dynamic_cast<const SMDS_VtkVolume*>( elem );
-            if (aPolyedre)
-            {
-              int nbFaces = aPolyedre->NbFaces();
-
-              vector<const SMDS_MeshNode *> poly_nodes;
-              vector<int> quantities;
-
-              for (int iface = 1; iface <= nbFaces; iface++) {
-                int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
-                vector<const SMDS_MeshNode *> faceNodes (nbFaceNodes);
-
-                for (int inode = 1; inode <= nbFaceNodes; inode++) {
-                  const SMDS_MeshNode * faceNode = aPolyedre->GetFaceNode(iface, inode);
-                  TNodeNodeMap::iterator nnIt = nodeNodeMap.find(faceNode);
-                  if (nnIt != nodeNodeMap.end()) { // faceNode sticks
-                    faceNode = (*nnIt).second;
-                  }
-                  faceNodes[inode - 1] = faceNode;
-                }
+          bool isValid = ( nbNewNodes % 2 == 0 );
+          for ( int i = 0; i < nbNewNodes && isValid; ++i )
+            isValid = ( elem->IsMediumNode( face_nodes[i]) == bool( i % 2 ));
+          elemType->SetQuad( isValid );
+          if ( isValid ) // put medium nodes after corners
+            SMDS_MeshCell::applyInterlaceRev
+              ( SMDS_MeshCell::interlacedSmdsOrder( SMDSEntity_Quad_Polygon,
+                                                    nbNewNodes ), face_nodes );
+        }
+        elemType->SetPoly(( nbNewNodes / ( elemType->myIsQuad + 1 ) > 4 ));
+      }
+      nbUniqueNodes = newElemDefs[0].myNodes.size();
+      break;
+    } // Polygon
 
-                SimplifyFace(faceNodes, poly_nodes, quantities);
-              }
+    case SMDSEntity_Polyhedra: // Polyhedral volume
+    {
+      if ( nbUniqueNodes >= 4 )
+      {
+        // each face has to be analyzed in order to check volume validity
+        if ( const SMDS_MeshVolume* aPolyedre = SMDS_Mesh::DownCast< SMDS_MeshVolume >( elem ))
+        {
+          int nbFaces = aPolyedre->NbFaces();
 
-              if (quantities.size() > 3) {
-                // to be done: remove coincident faces
-              }
+          vector<const SMDS_MeshNode *>& poly_nodes = newElemDefs[0].myNodes;
+          vector<int>                  & quantities = newElemDefs[0].myPolyhedQuantities;
+          vector<const SMDS_MeshNode *>  faceNodes;
+          poly_nodes.clear();
+          quantities.clear();
 
-              if (quantities.size() > 3)
-              {
-                const SMDS_MeshElement* newElem =
-                  aMesh->AddPolyhedralVolume(poly_nodes, quantities);
-                myLastCreatedElems.Append(newElem);
-                if ( aShapeId && newElem )
-                  aMesh->SetMeshElementOnShape( newElem, aShapeId );
-                rmElemIds.push_back(elem->GetID());
-              }
-            }
-            else {
-              rmElemIds.push_back(elem->GetID());
+          for (int iface = 1; iface <= nbFaces; iface++)
+          {
+            int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
+            faceNodes.resize( nbFaceNodes );
+            for (int inode = 1; inode <= nbFaceNodes; inode++)
+            {
+              const SMDS_MeshNode * faceNode = aPolyedre->GetFaceNode(iface, inode);
+              TNodeNodeMap::iterator nnIt = nodeNodeMap.find(faceNode);
+              if ( nnIt != nodeNodeMap.end() ) // faceNode sticks
+                faceNode = (*nnIt).second;
+              faceNodes[inode - 1] = faceNode;
             }
+            SimplifyFace(faceNodes, poly_nodes, quantities);
           }
-        }
-        else {
-        }
 
-        continue;
-      } // poly element
-
-      // Regular elements
-      // TODO not all the possible cases are solved. Find something more generic?
-      switch ( nbNodes ) {
-      case 2: ///////////////////////////////////// EDGE
-        isOk = false; break;
-      case 3: ///////////////////////////////////// TRIANGLE
-        isOk = false; break;
-      case 4:
-        if ( elem->GetType() == SMDSAbs_Volume ) // TETRAHEDRON
-          isOk = false;
-        else { //////////////////////////////////// QUADRANGLE
-          if ( nbUniqueNodes < 3 )
-            isOk = false;
-          else if ( nbRepl == 2 && iRepl[ 1 ] - iRepl[ 0 ] == 2 )
-            isOk = false; // opposite nodes stick
-          //MESSAGE("isOk " << isOk);
-        }
-        break;
-      case 6: ///////////////////////////////////// PENTAHEDRON
-        if ( nbUniqueNodes == 4 ) {
-          // ---------------------------------> tetrahedron
-          if (nbRepl == 3 &&
-              iRepl[ 0 ] > 2 && iRepl[ 1 ] > 2 && iRepl[ 2 ] > 2 ) {
-            // all top nodes stick: reverse a bottom
-            uniqueNodes[ 0 ] = curNodes [ 1 ];
-            uniqueNodes[ 1 ] = curNodes [ 0 ];
-          }
-          else if (nbRepl == 3 &&
-                   iRepl[ 0 ] < 3 && iRepl[ 1 ] < 3 && iRepl[ 2 ] < 3 ) {
-            // all bottom nodes stick: set a top before
-            uniqueNodes[ 3 ] = uniqueNodes [ 0 ];
-            uniqueNodes[ 0 ] = curNodes [ 3 ];
-            uniqueNodes[ 1 ] = curNodes [ 4 ];
-            uniqueNodes[ 2 ] = curNodes [ 5 ];
-          }
-          else if (nbRepl == 4 &&
-                   iRepl[ 2 ] - iRepl [ 0 ] == 3 && iRepl[ 3 ] - iRepl [ 1 ] == 3 ) {
-            // a lateral face turns into a line: reverse a bottom
-            uniqueNodes[ 0 ] = curNodes [ 1 ];
-            uniqueNodes[ 1 ] = curNodes [ 0 ];
-          }
-          else
-            isOk = false;
-        }
-        else if ( nbUniqueNodes == 5 ) {
-          // PENTAHEDRON --------------------> 2 tetrahedrons
-          if ( nbRepl == 2 && iRepl[ 1 ] - iRepl [ 0 ] == 3 ) {
-            // a bottom node sticks with a linked top one
-            // 1.
-            SMDS_MeshElement* newElem =
-              aMesh->AddVolume(curNodes[ 3 ],
-                               curNodes[ 4 ],
-                               curNodes[ 5 ],
-                               curNodes[ iRepl[ 0 ] == 2 ? 1 : 2 ]);
-            myLastCreatedElems.Append(newElem);
-            if ( aShapeId )
-              aMesh->SetMeshElementOnShape( newElem, aShapeId );
-            // 2. : reverse a bottom
-            uniqueNodes[ 0 ] = curNodes [ 1 ];
-            uniqueNodes[ 1 ] = curNodes [ 0 ];
-            nbUniqueNodes = 4;
+          if ( quantities.size() > 3 )
+          {
+            // TODO: remove coincident faces
+            nbResElems = 1;
+            nbUniqueNodes = newElemDefs[0].myNodes.size();
           }
-          else
-            isOk = false;
         }
-        else
-          isOk = false;
-        break;
-      case 8: {
-        if(elem->IsQuadratic()) { // Quadratic quadrangle
-          //   1    5    2
-          //    +---+---+
-          //    |       |
-          //    |       |
-          //   4+       +6
-          //    |       |
-          //    |       |
-          //    +---+---+
-          //   0    7    3
-          isOk = false;
-          if(nbRepl==2) {
-            MESSAGE("nbRepl=2: " << iRepl[0] << " " << iRepl[1]);
-          }
-          if(nbRepl==3) {
-            MESSAGE("nbRepl=3: " << iRepl[0] << " " << iRepl[1]  << " " << iRepl[2]);
-            nbUniqueNodes = 6;
-            if( iRepl[0]==0 && iRepl[1]==1 && iRepl[2]==4 ) {
-              uniqueNodes[0] = curNodes[0];
-              uniqueNodes[1] = curNodes[2];
-              uniqueNodes[2] = curNodes[3];
-              uniqueNodes[3] = curNodes[5];
-              uniqueNodes[4] = curNodes[6];
-              uniqueNodes[5] = curNodes[7];
-              isOk = true;
-            }
-            if( iRepl[0]==0 && iRepl[1]==3 && iRepl[2]==7 ) {
-              uniqueNodes[0] = curNodes[0];
-              uniqueNodes[1] = curNodes[1];
-              uniqueNodes[2] = curNodes[2];
-              uniqueNodes[3] = curNodes[4];
-              uniqueNodes[4] = curNodes[5];
-              uniqueNodes[5] = curNodes[6];
-              isOk = true;
-            }
-            if( iRepl[0]==0 && iRepl[1]==4 && iRepl[2]==7 ) {
-              uniqueNodes[0] = curNodes[1];
-              uniqueNodes[1] = curNodes[2];
-              uniqueNodes[2] = curNodes[3];
-              uniqueNodes[3] = curNodes[5];
-              uniqueNodes[4] = curNodes[6];
-              uniqueNodes[5] = curNodes[0];
-              isOk = true;
-            }
-            if( iRepl[0]==1 && iRepl[1]==2 && iRepl[2]==5 ) {
-              uniqueNodes[0] = curNodes[0];
-              uniqueNodes[1] = curNodes[1];
-              uniqueNodes[2] = curNodes[3];
-              uniqueNodes[3] = curNodes[4];
-              uniqueNodes[4] = curNodes[6];
-              uniqueNodes[5] = curNodes[7];
-              isOk = true;
-            }
-            if( iRepl[0]==1 && iRepl[1]==4 && iRepl[2]==5 ) {
-              uniqueNodes[0] = curNodes[0];
-              uniqueNodes[1] = curNodes[2];
-              uniqueNodes[2] = curNodes[3];
-              uniqueNodes[3] = curNodes[1];
-              uniqueNodes[4] = curNodes[6];
-              uniqueNodes[5] = curNodes[7];
-              isOk = true;
-            }
-            if( iRepl[0]==2 && iRepl[1]==3 && iRepl[2]==6 ) {
-              uniqueNodes[0] = curNodes[0];
-              uniqueNodes[1] = curNodes[1];
-              uniqueNodes[2] = curNodes[2];
-              uniqueNodes[3] = curNodes[4];
-              uniqueNodes[4] = curNodes[5];
-              uniqueNodes[5] = curNodes[7];
-              isOk = true;
-            }
-            if( iRepl[0]==2 && iRepl[1]==5 && iRepl[2]==6 ) {
-              uniqueNodes[0] = curNodes[0];
-              uniqueNodes[1] = curNodes[1];
-              uniqueNodes[2] = curNodes[3];
-              uniqueNodes[3] = curNodes[4];
-              uniqueNodes[4] = curNodes[2];
-              uniqueNodes[5] = curNodes[7];
-              isOk = true;
+      }
+    }
+    break;
+
+    // Regular elements
+    // TODO not all the possible cases are solved. Find something more generic?
+    case SMDSEntity_Edge: //////// EDGE
+    case SMDSEntity_Triangle: //// TRIANGLE
+    case SMDSEntity_Quad_Triangle:
+    case SMDSEntity_Tetra:
+    case SMDSEntity_Quad_Tetra: // TETRAHEDRON
+    {
+      break;
+    }
+    case SMDSEntity_Quad_Edge:
+    {
+      break;
+    }
+    case SMDSEntity_Quadrangle: //////////////////////////////////// QUADRANGLE
+    {
+      if ( nbUniqueNodes < 3 )
+        toRemove = true;
+      else if ( nbRepl == 1 && curNodes[ iRepl[0]] == curNodes[( iRepl[0]+2 )%4 ])
+        toRemove = true; // opposite nodes stick
+      else
+        toRemove = false;
+      break;
+    }
+    case SMDSEntity_Quad_Quadrangle: // Quadratic QUADRANGLE
+    {
+      //   1    5    2
+      //    +---+---+
+      //    |       |
+      //   4+       +6
+      //    |       |
+      //    +---+---+
+      //   0    7    3
+      if ( nbUniqueNodes == 6 &&
+           iRepl[0] < 4       &&
+           ( nbRepl == 1 || iRepl[1] >= 4 ))
+      {
+        toRemove = false;
+      }
+      break;
+    }
+    case SMDSEntity_BiQuad_Quadrangle: // Bi-Quadratic QUADRANGLE
+    {
+      //   1    5    2
+      //    +---+---+
+      //    |       |
+      //   4+  8+   +6
+      //    |       |
+      //    +---+---+
+      //   0    7    3
+      if ( nbUniqueNodes == 7 &&
+           iRepl[0] < 4       &&
+           ( nbRepl == 1 || iRepl[1] != 8 ))
+      {
+        toRemove = false;
+      }
+      break;
+    }
+    case SMDSEntity_Penta: ///////////////////////////////////// PENTAHEDRON
+    {
+      if ( nbUniqueNodes == 4 ) {
+        // ---------------------------------> tetrahedron
+        if ( curNodes[3] == curNodes[4] &&
+             curNodes[3] == curNodes[5] ) {
+          // top nodes stick
+          toRemove = false;
+        }
+        else if ( curNodes[0] == curNodes[1] &&
+                  curNodes[0] == curNodes[2] ) {
+          // bottom nodes stick: set a top before
+          uniqueNodes[ 3 ] = uniqueNodes [ 0 ];
+          uniqueNodes[ 0 ] = curNodes [ 5 ];
+          uniqueNodes[ 1 ] = curNodes [ 4 ];
+          uniqueNodes[ 2 ] = curNodes [ 3 ];
+          toRemove = false;
+        }
+        else if (( curNodes[0] == curNodes[3] ) +
+                 ( curNodes[1] == curNodes[4] ) +
+                 ( curNodes[2] == curNodes[5] ) == 2 ) {
+          // a lateral face turns into a line
+          toRemove = false;
+        }
+      }
+      else if ( nbUniqueNodes == 5 ) {
+        // PENTAHEDRON --------------------> pyramid
+        if ( curNodes[0] == curNodes[3] )
+        {
+          uniqueNodes[ 0 ] = curNodes[ 1 ];
+          uniqueNodes[ 1 ] = curNodes[ 4 ];
+          uniqueNodes[ 2 ] = curNodes[ 5 ];
+          uniqueNodes[ 3 ] = curNodes[ 2 ];
+          uniqueNodes[ 4 ] = curNodes[ 0 ];
+          toRemove = false;
+        }
+        if ( curNodes[1] == curNodes[4] )
+        {
+          uniqueNodes[ 0 ] = curNodes[ 0 ];
+          uniqueNodes[ 1 ] = curNodes[ 2 ];
+          uniqueNodes[ 2 ] = curNodes[ 5 ];
+          uniqueNodes[ 3 ] = curNodes[ 3 ];
+          uniqueNodes[ 4 ] = curNodes[ 1 ];
+          toRemove = false;
+        }
+        if ( curNodes[2] == curNodes[5] )
+        {
+          uniqueNodes[ 0 ] = curNodes[ 0 ];
+          uniqueNodes[ 1 ] = curNodes[ 3 ];
+          uniqueNodes[ 2 ] = curNodes[ 4 ];
+          uniqueNodes[ 3 ] = curNodes[ 1 ];
+          uniqueNodes[ 4 ] = curNodes[ 2 ];
+          toRemove = false;
+        }
+      }
+      break;
+    }
+    case SMDSEntity_Hexa:
+    {
+      //////////////////////////////////// HEXAHEDRON
+      SMDS_VolumeTool hexa (elem);
+      hexa.SetExternalNormal();
+      if ( nbUniqueNodes == 4 && nbRepl == 4 ) {
+        //////////////////////// HEX ---> tetrahedron
+        for ( int iFace = 0; iFace < 6; iFace++ ) {
+          const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
+          if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
+              curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
+              curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
+            // one face turns into a point ...
+            int  pickInd = ind[ 0 ];
+            int iOppFace = hexa.GetOppFaceIndex( iFace );
+            ind = hexa.GetFaceNodesIndices( iOppFace );
+            int nbStick = 0;
+            uniqueNodes.clear();
+            for ( iCur = 0; iCur < 4 && nbStick < 2; iCur++ ) {
+              if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
+                nbStick++;
+              else
+                uniqueNodes.push_back( curNodes[ind[ iCur ]]);
             }
-            if( iRepl[0]==3 && iRepl[1]==6 && iRepl[2]==7 ) {
-              uniqueNodes[0] = curNodes[0];
-              uniqueNodes[1] = curNodes[1];
-              uniqueNodes[2] = curNodes[2];
-              uniqueNodes[3] = curNodes[4];
-              uniqueNodes[4] = curNodes[5];
-              uniqueNodes[5] = curNodes[3];
-              isOk = true;
+            if ( nbStick == 1 ) {
+              // ... and the opposite one - into a triangle.
+              // set a top node
+              uniqueNodes.push_back( curNodes[ pickInd ]);
+              toRemove = false;
             }
+            break;
           }
-          if(nbRepl==4) {
-            MESSAGE("nbRepl=4: " << iRepl[0] << " " << iRepl[1]  << " " << iRepl[2] << " " << iRepl[3]);
-          }
-          if(nbRepl==5) {
-            MESSAGE("nbRepl=5: " << iRepl[0] << " " << iRepl[1]  << " " << iRepl[2] << " " << iRepl[3] << " " << iRepl[4]);
-          }
-          break;
         }
-        //////////////////////////////////// HEXAHEDRON
-        isOk = false;
-        SMDS_VolumeTool hexa (elem);
-        hexa.SetExternalNormal();
-        if ( nbUniqueNodes == 4 && nbRepl == 4 ) {
-          //////////////////////// HEX ---> 1 tetrahedron
-          for ( int iFace = 0; iFace < 6; iFace++ ) {
-            const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
-            if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
-                curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
-                curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
-              // one face turns into a point ...
-              int iOppFace = hexa.GetOppFaceIndex( iFace );
-              ind = hexa.GetFaceNodesIndices( iOppFace );
-              int nbStick = 0;
-              for ( iCur = 0; iCur < 4 && nbStick < 2; iCur++ ) {
-                if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
-                  nbStick++;
-              }
-              if ( nbStick == 1 ) {
-                // ... and the opposite one - into a triangle.
-                // set a top node
-                ind = hexa.GetFaceNodesIndices( iFace );
-                uniqueNodes[ 3 ] = curNodes[ind[ 0 ]];
-                isOk = true;
+      }
+      else if ( nbUniqueNodes == 6 && nbRepl == 2 ) {
+        //////////////////////// HEX ---> prism
+        int nbTria = 0, iTria[3];
+        const int *ind; // indices of face nodes
+        // look for triangular faces
+        for ( int iFace = 0; iFace < 6 && nbTria < 3; iFace++ ) {
+          ind = hexa.GetFaceNodesIndices( iFace );
+          TIDSortedNodeSet faceNodes;
+          for ( iCur = 0; iCur < 4; iCur++ )
+            faceNodes.insert( curNodes[ind[iCur]] );
+          if ( faceNodes.size() == 3 )
+            iTria[ nbTria++ ] = iFace;
+        }
+        // check if triangles are opposite
+        if ( nbTria == 2 && iTria[0] == hexa.GetOppFaceIndex( iTria[1] ))
+        {
+          // set nodes of the bottom triangle
+          ind = hexa.GetFaceNodesIndices( iTria[ 0 ]);
+          vector<int> indB;
+          for ( iCur = 0; iCur < 4; iCur++ )
+            if ( ind[iCur] != iRepl[0] && ind[iCur] != iRepl[1])
+              indB.push_back( ind[iCur] );
+          if ( !hexa.IsForward() )
+            std::swap( indB[0], indB[2] );
+          for ( iCur = 0; iCur < 3; iCur++ )
+            uniqueNodes[ iCur ] = curNodes[indB[iCur]];
+          // set nodes of the top triangle
+          const int *indT = hexa.GetFaceNodesIndices( iTria[ 1 ]);
+          for ( iCur = 0; iCur < 3; ++iCur )
+            for ( int j = 0; j < 4; ++j )
+              if ( hexa.IsLinked( indB[ iCur ], indT[ j ] ))
+              {
+                uniqueNodes[ iCur + 3 ] = curNodes[ indT[ j ]];
+                break;
               }
-              break;
-            }
-          }
-        }
-        else if ( nbUniqueNodes == 6 && nbRepl == 2 ) {
-          //////////////////////// HEX ---> 1 prism
-          int nbTria = 0, iTria[3];
-          const int *ind; // indices of face nodes
-          // look for triangular faces
-          for ( int iFace = 0; iFace < 6 && nbTria < 3; iFace++ ) {
-            ind = hexa.GetFaceNodesIndices( iFace );
-            TIDSortedNodeSet faceNodes;
-            for ( iCur = 0; iCur < 4; iCur++ )
-              faceNodes.insert( curNodes[ind[iCur]] );
-            if ( faceNodes.size() == 3 )
-              iTria[ nbTria++ ] = iFace;
-          }
-          // check if triangles are opposite
-          if ( nbTria == 2 && iTria[0] == hexa.GetOppFaceIndex( iTria[1] ))
-          {
-            isOk = true;
-            // set nodes of the bottom triangle
-            ind = hexa.GetFaceNodesIndices( iTria[ 0 ]);
-            vector<int> indB;
-            for ( iCur = 0; iCur < 4; iCur++ )
-              if ( ind[iCur] != iRepl[0] && ind[iCur] != iRepl[1])
-                indB.push_back( ind[iCur] );
-            if ( !hexa.IsForward() )
-              std::swap( indB[0], indB[2] );
-            for ( iCur = 0; iCur < 3; iCur++ )
-              uniqueNodes[ iCur ] = curNodes[indB[iCur]];
-            // set nodes of the top triangle
-            const int *indT = hexa.GetFaceNodesIndices( iTria[ 1 ]);
-            for ( iCur = 0; iCur < 3; ++iCur )
-              for ( int j = 0; j < 4; ++j )
-                if ( hexa.IsLinked( indB[ iCur ], indT[ j ] ))
-                {
-                  uniqueNodes[ iCur + 3 ] = curNodes[ indT[ j ]];
-                  break;
-                }
-          }
+          toRemove = false;
           break;
         }
-        else if (nbUniqueNodes == 5 && nbRepl == 4 ) {
-          //////////////////// HEXAHEDRON ---> 2 tetrahedrons
-          for ( int iFace = 0; iFace < 6; iFace++ ) {
-            const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
-            if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
-                curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
-                curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
-              // one face turns into a point ...
-              int iOppFace = hexa.GetOppFaceIndex( iFace );
-              ind = hexa.GetFaceNodesIndices( iOppFace );
-              int nbStick = 0;
-              iUnique = 2;  // reverse a tetrahedron 1 bottom
-              for ( iCur = 0; iCur < 4 && nbStick == 0; iCur++ ) {
-                if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
-                  nbStick++;
-                else if ( iUnique >= 0 )
-                  uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
-              }
-              if ( nbStick == 0 ) {
-                // ... and the opposite one is a quadrangle
-                // set a top node
-                const int* indTop = hexa.GetFaceNodesIndices( iFace );
-                uniqueNodes[ 3 ] = curNodes[indTop[ 0 ]];
-                nbUniqueNodes = 4;
-                // tetrahedron 2
-                SMDS_MeshElement* newElem =
-                  aMesh->AddVolume(curNodes[ind[ 0 ]],
-                                   curNodes[ind[ 3 ]],
-                                   curNodes[ind[ 2 ]],
-                                   curNodes[indTop[ 0 ]]);
-                myLastCreatedElems.Append(newElem);
-                if ( aShapeId )
-                  aMesh->SetMeshElementOnShape( newElem, aShapeId );
-                isOk = true;
-              }
-              break;
+      }
+      else if (nbUniqueNodes == 5 && nbRepl == 3 ) {
+        //////////////////// HEXAHEDRON ---> pyramid
+        for ( int iFace = 0; iFace < 6; iFace++ ) {
+          const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
+          if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
+              curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
+              curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
+            // one face turns into a point ...
+            int iOppFace = hexa.GetOppFaceIndex( iFace );
+            ind = hexa.GetFaceNodesIndices( iOppFace );
+            uniqueNodes.clear();
+            for ( iCur = 0; iCur < 4; iCur++ ) {
+              if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
+                break;
+              else
+                uniqueNodes.push_back( curNodes[ind[ iCur ]]);
             }
+            if ( uniqueNodes.size() == 4 ) {
+              // ... and the opposite one is a quadrangle
+              // set a top node
+              const int* indTop = hexa.GetFaceNodesIndices( iFace );
+              uniqueNodes.push_back( curNodes[indTop[ 0 ]]);
+              toRemove = false;
+            }
+            break;
           }
         }
-        else if ( nbUniqueNodes == 6 && nbRepl == 4 ) {
-          ////////////////// HEXAHEDRON ---> 2 tetrahedrons or 1 prism
-          // find indices of quad and tri faces
-          int iQuadFace[ 6 ], iTriFace[ 6 ], nbQuad = 0, nbTri = 0, iFace;
-          for ( iFace = 0; iFace < 6; iFace++ ) {
-            const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
-            nodeSet.clear();
-            for ( iCur = 0; iCur < 4; iCur++ )
-              nodeSet.insert( curNodes[ind[ iCur ]] );
-            nbUniqueNodes = nodeSet.size();
-            if ( nbUniqueNodes == 3 )
-              iTriFace[ nbTri++ ] = iFace;
-            else if ( nbUniqueNodes == 4 )
-              iQuadFace[ nbQuad++ ] = iFace;
-          }
-          if (nbQuad == 2 && nbTri == 4 &&
-              hexa.GetOppFaceIndex( iQuadFace[ 0 ] ) == iQuadFace[ 1 ]) {
-            // 2 opposite quadrangles stuck with a diagonal;
-            // sample groups of merged indices: (0-4)(2-6)
-            // --------------------------------------------> 2 tetrahedrons
-            const int *ind1 = hexa.GetFaceNodesIndices( iQuadFace[ 0 ]); // indices of quad1 nodes
-            const int *ind2 = hexa.GetFaceNodesIndices( iQuadFace[ 1 ]);
-            int i0, i1d, i2, i3d, i0t, i2t; // d-daigonal, t-top
-            if (curNodes[ind1[ 0 ]] == curNodes[ind2[ 0 ]] &&
-                curNodes[ind1[ 2 ]] == curNodes[ind2[ 2 ]]) {
-              // stuck with 0-2 diagonal
-              i0  = ind1[ 3 ];
-              i1d = ind1[ 0 ];
-              i2  = ind1[ 1 ];
-              i3d = ind1[ 2 ];
-              i0t = ind2[ 1 ];
-              i2t = ind2[ 3 ];
-            }
-            else if (curNodes[ind1[ 1 ]] == curNodes[ind2[ 3 ]] &&
-                     curNodes[ind1[ 3 ]] == curNodes[ind2[ 1 ]]) {
-              // stuck with 1-3 diagonal
-              i0  = ind1[ 0 ];
-              i1d = ind1[ 1 ];
-              i2  = ind1[ 2 ];
-              i3d = ind1[ 3 ];
-              i0t = ind2[ 0 ];
-              i2t = ind2[ 1 ];
-            }
-            else {
-              ASSERT(0);
-            }
-            // tetrahedron 1
-            uniqueNodes[ 0 ] = curNodes [ i0 ];
-            uniqueNodes[ 1 ] = curNodes [ i1d ];
-            uniqueNodes[ 2 ] = curNodes [ i3d ];
-            uniqueNodes[ 3 ] = curNodes [ i0t ];
-            nbUniqueNodes = 4;
-            // tetrahedron 2
-            SMDS_MeshElement* newElem = aMesh->AddVolume(curNodes[ i1d ],
-                                                         curNodes[ i2 ],
-                                                         curNodes[ i3d ],
-                                                         curNodes[ i2t ]);
-            myLastCreatedElems.Append(newElem);
-            if ( aShapeId )
-              aMesh->SetMeshElementOnShape( newElem, aShapeId );
-            isOk = true;
+      }
+
+      if ( toRemove && nbUniqueNodes > 4 ) {
+        ////////////////// HEXAHEDRON ---> polyhedron
+        hexa.SetExternalNormal();
+        vector<const SMDS_MeshNode *>& poly_nodes = newElemDefs[0].myNodes;
+        vector<int>                  & quantities = newElemDefs[0].myPolyhedQuantities;
+        poly_nodes.reserve( 6 * 4 ); poly_nodes.clear();
+        quantities.reserve( 6 );     quantities.clear();
+        for ( int iFace = 0; iFace < 6; iFace++ )
+        {
+          const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
+          if ( curNodes[ind[0]] == curNodes[ind[2]] ||
+               curNodes[ind[1]] == curNodes[ind[3]] )
+          {
+            quantities.clear();
+            break; // opposite nodes stick
           }
-          else if (( nbTri == 2 && nbQuad == 3 ) || // merged (0-4)(1-5)
-                   ( nbTri == 4 && nbQuad == 2 )) { // merged (7-4)(1-5)
-            // --------------------------------------------> prism
-            // find 2 opposite triangles
-            nbUniqueNodes = 6;
-            for ( iFace = 0; iFace + 1 < nbTri; iFace++ ) {
-              if ( hexa.GetOppFaceIndex( iTriFace[ iFace ] ) == iTriFace[ iFace + 1 ]) {
-                // find indices of kept and replaced nodes
-                // and fill unique nodes of 2 opposite triangles
-                const int *ind1 = hexa.GetFaceNodesIndices( iTriFace[ iFace ]);
-                const int *ind2 = hexa.GetFaceNodesIndices( iTriFace[ iFace + 1 ]);
-                const SMDS_MeshNode** hexanodes = hexa.GetNodes();
-                // fill unique nodes
-                iUnique = 0;
-                isOk = true;
-                for ( iCur = 0; iCur < 4 && isOk; iCur++ ) {
-                  const SMDS_MeshNode* n     = curNodes[ind1[ iCur ]];
-                  const SMDS_MeshNode* nInit = hexanodes[ind1[ iCur ]];
-                  if ( n == nInit ) {
-                    // iCur of a linked node of the opposite face (make normals co-directed):
-                    int iCurOpp = ( iCur == 1 || iCur == 3 ) ? 4 - iCur : iCur;
-                    // check that correspondent corners of triangles are linked
-                    if ( !hexa.IsLinked( ind1[ iCur ], ind2[ iCurOpp ] ))
-                      isOk = false;
-                    else {
-                      uniqueNodes[ iUnique ] = n;
-                      uniqueNodes[ iUnique + 3 ] = curNodes[ind2[ iCurOpp ]];
-                      iUnique++;
-                    }
-                  }
-                }
-                break;
-              }
-            }
+          nodeSet.clear();
+          for ( iCur = 0; iCur < 4; iCur++ )
+          {
+            if ( nodeSet.insert( curNodes[ind[ iCur ]] ).second )
+              poly_nodes.push_back( curNodes[ind[ iCur ]]);
           }
-        } // if ( nbUniqueNodes == 6 && nbRepl == 4 )
-        else
+          if ( nodeSet.size() < 3 )
+            poly_nodes.resize( poly_nodes.size() - nodeSet.size() );
+          else
+            quantities.push_back( nodeSet.size() );
+        }
+        if ( quantities.size() >= 4 )
         {
-          MESSAGE("MergeNodes() removes hexahedron "<< elem);
+          nbResElems = 1;
+          nbUniqueNodes = poly_nodes.size();
+          newElemDefs[0].SetPoly(true);
         }
-        break;
-      } // HEXAHEDRON
+      }
+      break;
+    } // case HEXAHEDRON
 
-      default:
-        isOk = false;
-      } // switch ( nbNodes )
+    default:
+      toRemove = true;
 
-    } // if ( nbNodes != nbUniqueNodes ) // some nodes stick
+    } // switch ( entity )
 
-    if ( isOk ) // the non-poly elem remains valid after sticking nodes
+    if ( toRemove && nbResElems == 0 && avoidMakingHoles )
     {
-      if ( nbNodes != nbUniqueNodes ||
-           !aMesh->ChangeElementNodes( elem, & curNodes[0], nbNodes ))
-      {
-        elemType.Init( elem ).SetID( elem->GetID() );
-
-        SMESHDS_SubMesh * sm = aShapeId > 0 ? aMesh->MeshElements(aShapeId) : 0;
-        aMesh->RemoveFreeElement(elem, sm, /*fromGroups=*/false);
-
-        uniqueNodes.resize(nbUniqueNodes);
-        SMDS_MeshElement* newElem = this->AddElement( uniqueNodes, elemType );
-        if ( sm && newElem )
-          sm->AddElement( newElem );
-        if ( elem != newElem )
-          ReplaceElemInGroups( elem, newElem, aMesh );
-      }
-    }
-    else {
-      // Remove invalid regular element or invalid polygon
-      rmElemIds.push_back( elem->GetID() );
+      // erase from nodeNodeMap nodes whose merge spoils elem
+      vector< const SMDS_MeshNode* > noMergeNodes;
+      SMESH_MeshAlgos::DeMerge( elem, curNodes, noMergeNodes );
+      for ( size_t i = 0; i < noMergeNodes.size(); ++i )
+        nodeNodeMap.erase( noMergeNodes[i] );
     }
+    
+  } // if ( nbNodes != nbUniqueNodes ) // some nodes stick
 
-  } // loop on elements
+  uniqueNodes.resize( nbUniqueNodes );
 
-  // Remove bad elements, then equal nodes (order important)
+  if ( !toRemove && nbResElems == 0 )
+    nbResElems = 1;
 
-  Remove( rmElemIds, false );
-  Remove( rmNodeIds, true );
+  newElemDefs.resize( nbResElems );
 
-  return;
+  return !toRemove;
 }
 
 
 // ========================================================
-// class   : SortableElement
-// purpose : allow sorting elements basing on their nodes
+// class   : ComparableElement
+// purpose : allow comparing elements basing on their nodes
 // ========================================================
-class SortableElement : public set <const SMDS_MeshElement*>
+
+class ComparableElement : public boost::container::flat_set< int >
 {
+  typedef boost::container::flat_set< int >  int_set;
+
+  const SMDS_MeshElement* myElem;
+  int                     mySumID;
+  mutable int             myGroupID;
+
 public:
 
-  SortableElement( const SMDS_MeshElement* theElem )
+  ComparableElement( const SMDS_MeshElement* theElem ):
+    myElem ( theElem ), mySumID( 0 ), myGroupID( -1 )
   {
-    myElem = theElem;
-    SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
-    while ( nodeIt->more() )
-      this->insert( nodeIt->next() );
+    this->reserve( theElem->NbNodes() );
+    for ( SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator(); nodeIt->more(); )
+    {
+      int id = nodeIt->next()->GetID();
+      mySumID += id;
+      this->insert( id );
+    }
   }
 
-  const SMDS_MeshElement* Get() const
-  { return myElem; }
+  const SMDS_MeshElement* GetElem() const { return myElem; }
+
+  int& GroupID() const { return myGroupID; }
+  //int& GroupID() const { return const_cast< int& >( myGroupID ); }
+
+  ComparableElement( const ComparableElement& theSource ) // move copy
+  {
+    ComparableElement& src = const_cast< ComparableElement& >( theSource );
+    (int_set&) (*this ) = boost::move( src );
+    myElem    = src.myElem;
+    mySumID   = src.mySumID;
+    myGroupID = src.myGroupID;
+  }
+
+  static int HashCode(const ComparableElement& se, int limit )
+  {
+    return ::HashCode( se.mySumID, limit );
+  }
+  static Standard_Boolean IsEqual(const ComparableElement& se1, const ComparableElement& se2 )
+  {
+    return ( se1 == se2 );
+  }
 
-private:
-  mutable const SMDS_MeshElement* myElem;
 };
 
 //=======================================================================
@@ -7966,53 +7799,47 @@ private:
 //           Search among theElements or in the whole mesh if theElements is empty
 //=======================================================================
 
-void SMESH_MeshEditor::FindEqualElements(TIDSortedElemSet &        theElements,
-                                         TListOfListOfElementsID & theGroupsOfElementsID)
+void SMESH_MeshEditor::FindEqualElements( TIDSortedElemSet &        theElements,
+                                          TListOfListOfElementsID & theGroupsOfElementsID )
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
-
-  typedef map< SortableElement, int > TMapOfNodeSet;
-  typedef list<int> TGroupOfElems;
+  ClearLastCreated();
 
-  if ( theElements.empty() )
-  { // get all elements in the mesh
-    SMDS_ElemIteratorPtr eIt = GetMeshDS()->elementsIterator();
-    while ( eIt->more() )
-      theElements.insert( theElements.end(), eIt->next() );
-  }
+  SMDS_ElemIteratorPtr elemIt;
+  if ( theElements.empty() ) elemIt = GetMeshDS()->elementsIterator();
+  else                       elemIt = SMESHUtils::elemSetIterator( theElements );
 
-  vector< TGroupOfElems > arrayOfGroups;
-  TGroupOfElems groupOfElems;
-  TMapOfNodeSet mapOfNodeSet;
+  typedef NCollection_Map< ComparableElement, ComparableElement > TMapOfElements;
+  typedef std::list<int>                                          TGroupOfElems;
+  TMapOfElements               mapOfElements;
+  std::vector< TGroupOfElems > arrayOfGroups;
+  TGroupOfElems                groupOfElems;
 
-  TIDSortedElemSet::iterator elemIt = theElements.begin();
-  for ( int i = 0; elemIt != theElements.end(); ++elemIt )
+  while ( elemIt->more() )
   {
-    const SMDS_MeshElement* curElem = *elemIt;
-    SortableElement SE(curElem);
+    const SMDS_MeshElement* curElem = elemIt->next();
+    ComparableElement      compElem = curElem;
     // check uniqueness
-    pair< TMapOfNodeSet::iterator, bool> pp = mapOfNodeSet.insert(make_pair(SE, i));
-    if ( !pp.second ) { // one more coincident elem
-      TMapOfNodeSet::iterator& itSE = pp.first;
-      int ind = (*itSE).second;
-      arrayOfGroups[ind].push_back( curElem->GetID() );
-    }
-    else {
-      arrayOfGroups.push_back( groupOfElems );
-      arrayOfGroups.back().push_back( curElem->GetID() );
-      i++;
+    const ComparableElement& elemInSet = mapOfElements.Added( compElem );
+    if ( elemInSet.GetElem() != curElem ) // coincident elem
+    {
+      int& iG = elemInSet.GroupID();
+      if ( iG < 0 )
+      {
+        iG = arrayOfGroups.size();
+        arrayOfGroups.push_back( groupOfElems );
+        arrayOfGroups[ iG ].push_back( elemInSet.GetElem()->GetID() );
+      }
+      arrayOfGroups[ iG ].push_back( curElem->GetID() );
     }
   }
 
   groupOfElems.clear();
-  vector< TGroupOfElems >::iterator groupIt = arrayOfGroups.begin();
+  std::vector< TGroupOfElems >::iterator groupIt = arrayOfGroups.begin();
   for ( ; groupIt != arrayOfGroups.end(); ++groupIt )
   {
     if ( groupIt->size() > 1 ) {
-      //groupOfElems.sort(); -- theElements is sorted already
-      theGroupsOfElementsID.push_back( groupOfElems );
-      theGroupsOfElementsID.back().splice( theGroupsOfElementsID.back().end(), *groupIt );
+      //groupOfElems.sort(); -- theElements are sorted already
+      theGroupsOfElementsID.emplace_back( *groupIt );
     }
   }
 }
@@ -8024,8 +7851,7 @@ void SMESH_MeshEditor::FindEqualElements(TIDSortedElemSet &        theElements,
 
 void SMESH_MeshEditor::MergeElements(TListOfListOfElementsID & theGroupsOfElementsID)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
   typedef list<int> TListOfIDs;
   TListOfIDs rmElemIds; // IDs of elems to remove
@@ -8064,8 +7890,8 @@ void SMESH_MeshEditor::MergeEqualElements()
   TIDSortedElemSet aMeshElements; /* empty input ==
                                      to merge equal elements in the whole mesh */
   TListOfListOfElementsID aGroupsOfElementsID;
-  FindEqualElements(aMeshElements, aGroupsOfElementsID);
-  MergeElements(aGroupsOfElementsID);
+  FindEqualElements( aMeshElements, aGroupsOfElementsID );
+  MergeElements( aGroupsOfElementsID );
 }
 
 //=======================================================================
@@ -8143,24 +7969,11 @@ bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode*             theFirst
       if ( e == curElem || foundElems.insert( e ).second ) {
         // get nodes
         int iNode = 0, nbNodes = e->NbNodes();
-        vector<const SMDS_MeshNode*> nodes(nbNodes+1);
-
-        if ( e->IsQuadratic() ) {
-          const SMDS_VtkFace* F =
-            dynamic_cast<const SMDS_VtkFace*>(e);
-          if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
-          // use special nodes iterator
-          SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
-          while( anIter->more() ) {
-            nodes[ iNode++ ] = cast2Node(anIter->next());
-          }
-        }
-        else {
-          SMDS_ElemIteratorPtr nIt = e->nodesIterator();
-          while ( nIt->more() )
-            nodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
-        }
-        nodes[ iNode ] = nodes[ 0 ];
+        vector<const SMDS_MeshNode*> nodes( nbNodes+1 );
+        nodes.assign( SMDS_MeshElement::iterator( e->interlacedNodesIterator() ),
+                      SMDS_MeshElement::iterator() );
+        nodes.push_back( nodes[ 0 ]);
+
         // check 2 links
         for ( iNode = 0; iNode < nbNodes; iNode++ )
           if (((nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) ||
@@ -8227,7 +8040,7 @@ bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode*             theFirst
       if ( contNodes[0].empty() && contNodes[1].empty() )
         return false;
 
-      // append the best free border
+      // push_back the best free border
       cNL = & contNodes[ contNodes[0].empty() ? 1 : 0 ];
       cFL = & contFaces[ contFaces[0].empty() ? 1 : 0 ];
       theNodes.pop_back(); // remove nIgnore
@@ -8277,10 +8090,8 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
                                  const bool           toCreatePolygons,
                                  const bool           toCreatePolyedrs)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
-  MESSAGE("::SewFreeBorder()");
   Sew_Error aResult = SEW_OK;
 
   // ====================================
@@ -8325,7 +8136,7 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
     // -------------------------------------------------------------------------
 
     // 1. Since sewing may break if there are volumes to split on the side 2,
-    //    we wont move nodes but just compute new coordinates for them
+    //    we won't move nodes but just compute new coordinates for them
     typedef map<const SMDS_MeshNode*, gp_XYZ> TNodeXYZMap;
     TNodeXYZMap nBordXYZ;
     list< const SMDS_MeshNode* >& bordNodes = nSide[ 0 ];
@@ -8388,7 +8199,7 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
     //const SMDS_MeshNode* faceNodes[ 4 ];
 
     const SMDS_MeshNode*    sideNode;
-    const SMDS_MeshElement* sideElem;
+    const SMDS_MeshElement* sideElem  = 0;
     const SMDS_MeshNode* prevSideNode = theSideFirstNode;
     const SMDS_MeshNode* prevBordNode = theBordFirstNode;
     nBordIt = bordNodes.begin();
@@ -8413,33 +8224,19 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
       {
         const SMDS_MeshElement* elem = invElemIt->next();
         // prepare data for a loop on links coming to prevSideNode, of a face or a volume
-        int iPrevNode, iNode = 0, nbNodes = elem->NbNodes();
+        int iPrevNode = 0, iNode = 0, nbNodes = elem->NbNodes();
         vector< const SMDS_MeshNode* > faceNodes( nbNodes, (const SMDS_MeshNode*)0 );
         bool isVolume = volume.Set( elem );
         const SMDS_MeshNode** nodes = isVolume ? volume.GetNodes() : & faceNodes[0];
         if ( isVolume ) // --volume
           hasVolumes = true;
-        else if ( elem->GetType()==SMDSAbs_Face ) { // --face
+        else if ( elem->GetType() == SMDSAbs_Face ) { // --face
           // retrieve all face nodes and find iPrevNode - an index of the prevSideNode
-          if(elem->IsQuadratic()) {
-            const SMDS_VtkFace* F =
-              dynamic_cast<const SMDS_VtkFace*>(elem);
-            if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
-            // use special nodes iterator
-            SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
-            while( anIter->more() ) {
-              nodes[ iNode ] = cast2Node(anIter->next());
-              if ( nodes[ iNode++ ] == prevSideNode )
-                iPrevNode = iNode - 1;
-            }
-          }
-          else {
-            SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
-            while ( nIt->more() ) {
-              nodes[ iNode ] = cast2Node( nIt->next() );
-              if ( nodes[ iNode++ ] == prevSideNode )
-                iPrevNode = iNode - 1;
-            }
+          SMDS_NodeIteratorPtr nIt = elem->interlacedNodesIterator();
+          while ( nIt->more() ) {
+            nodes[ iNode ] = cast2Node( nIt->next() );
+            if ( nodes[ iNode++ ] == prevSideNode )
+              iPrevNode = iNode - 1;
           }
           // there are 2 links to check
           nbNodes = 2;
@@ -8489,7 +8286,7 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
       } // loop on inverse elements of prevSideNode
 
       if ( !sideNode ) {
-        MESSAGE(" Cant find path by links of the Side 2 ");
+        MESSAGE(" Can't find path by links of the Side 2 ");
         return SEW_BAD_SIDE_NODES;
       }
       sideNodes.push_back( sideNode );
@@ -8732,19 +8529,20 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
   // get new segments
   TIDSortedElemSet segments;
   SMESH_SequenceOfElemPtr newFaces;
-  for ( int i = 1; i <= myLastCreatedElems.Length(); ++i )
+  for ( size_t i = 0; i < myLastCreatedElems.size(); ++i )
   {
-    if ( !myLastCreatedElems(i) ) continue;
-    if ( myLastCreatedElems(i)->GetType() == SMDSAbs_Edge )
-      segments.insert( segments.end(), myLastCreatedElems(i) );
+    if ( !myLastCreatedElems[i] ) continue;
+    if ( myLastCreatedElems[i]->GetType() == SMDSAbs_Edge )
+      segments.insert( segments.end(), myLastCreatedElems[i] );
     else
-      newFaces.Append( myLastCreatedElems(i) );
+      newFaces.push_back( myLastCreatedElems[i] );
   }
   // get segments adjacent to merged nodes
   TListOfListOfNodes::iterator groupIt = nodeGroupsToMerge.begin();
   for ( ; groupIt != nodeGroupsToMerge.end(); groupIt++ )
   {
     const list<const SMDS_MeshNode*>& nodes = *groupIt;
+    if ( nodes.front()->IsNull() ) continue;
     SMDS_ElemIteratorPtr segIt = nodes.front()->GetInverseElementIterator( SMDSAbs_Edge );
     while ( segIt->more() )
       segments.insert( segIt->next() );
@@ -8773,7 +8571,7 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
     myLastCreatedElems = newFaces;
     TIDSortedElemSet::iterator seg = segments.begin();
     for ( ; seg != segments.end(); ++seg )
-      myLastCreatedElems.Append( *seg );
+      myLastCreatedElems.push_back( *seg );
   }
 
   return aResult;
@@ -8843,7 +8641,7 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement*     theElemen
       if ( newElems[i] )
       {
         aMesh->SetMeshElementOnShape( newElems[i], theElement->getshapeId() );
-        myLastCreatedElems.Append( newElems[i] );
+        myLastCreatedElems.push_back( newElems[i] );
       }
     ReplaceElemInGroups( theElement, newElems, aMesh );
     aMesh->RemoveElement( theElement );
@@ -8899,48 +8697,21 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement*     theElemen
 
     // add nodes of face up to first node of link
     bool isFLN = false;
-
-    if ( theFace->IsQuadratic() ) {
-      const SMDS_VtkFace* F = dynamic_cast<const SMDS_VtkFace*>(theFace);
-      if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
-      // use special nodes iterator
-      SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
-      while( anIter->more()  && !isFLN ) {
-        const SMDS_MeshNode* n = cast2Node(anIter->next());
-        poly_nodes[iNode++] = n;
-        if (n == nodes[il1]) {
-          isFLN = true;
-        }
-      }
-      // add nodes to insert
-      list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
-      for (; nIt != aNodesToInsert.end(); nIt++) {
-        poly_nodes[iNode++] = *nIt;
-      }
-      // add nodes of face starting from last node of link
-      while ( anIter->more() ) {
-        poly_nodes[iNode++] = cast2Node(anIter->next());
-      }
+    SMDS_NodeIteratorPtr nodeIt = theFace->interlacedNodesIterator();
+    while ( nodeIt->more() && !isFLN ) {
+      const SMDS_MeshNode* n = nodeIt->next();
+      poly_nodes[iNode++] = n;
+      isFLN = ( n == nodes[il1] );
     }
-    else {
-      SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
-      while ( nodeIt->more() && !isFLN ) {
-        const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
-        poly_nodes[iNode++] = n;
-        if (n == nodes[il1]) {
-          isFLN = true;
-        }
-      }
-      // add nodes to insert
-      list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
-      for (; nIt != aNodesToInsert.end(); nIt++) {
-        poly_nodes[iNode++] = *nIt;
-      }
-      // add nodes of face starting from last node of link
-      while ( nodeIt->more() ) {
-        const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
-        poly_nodes[iNode++] = n;
-      }
+    // add nodes to insert
+    list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
+    for (; nIt != aNodesToInsert.end(); nIt++) {
+      poly_nodes[iNode++] = *nIt;
+    }
+    // add nodes of face starting from last node of link
+    while ( nodeIt->more() ) {
+      const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+      poly_nodes[iNode++] = n;
     }
 
     // make a new face
@@ -8961,7 +8732,7 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement*     theElemen
     }
     // decide how to split a quadrangle: compare possible variants
     // and choose which of splits to be a quadrangle
-    int i1, i2, iSplit, nbSplits = nbLinkNodes - 1, iBestQuad;
+    int i1, i2, iSplit, nbSplits = nbLinkNodes - 1, iBestQuad = 0;
     if ( nbFaceNodes == 3 ) {
       iBestQuad = nbSplits;
       i4 = i3;
@@ -9104,7 +8875,7 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement*     theElemen
     if ( newElems[i] )
     {
       aMesh->SetMeshElementOnShape( newElems[i], theFace->getshapeId() );
-      myLastCreatedElems.Append( newElems[i] );
+      myLastCreatedElems.push_back( newElems[i] );
     }
   ReplaceElemInGroups( theFace, newElems, aMesh );
   aMesh->RemoveElement(theFace);
@@ -9120,8 +8891,7 @@ void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode*        theBetweenNode
                                       const SMDS_MeshNode*        theBetweenNode2,
                                       list<const SMDS_MeshNode*>& theNodesToInsert)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
   SMDS_ElemIteratorPtr invElemIt = theBetweenNode1->GetInverseElementIterator(SMDSAbs_Volume);
   while (invElemIt->more()) { // loop on inverse elements of theBetweenNode1
@@ -9183,7 +8953,7 @@ void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode*        theBetweenNode
     if ( SMDS_MeshElement* newElem = aMesh->AddPolyhedralVolume( poly_nodes, quantities ))
     {
       aMesh->SetMeshElementOnShape( newElem, elem->getshapeId() );
-      myLastCreatedElems.Append( newElem );
+      myLastCreatedElems.push_back( newElem );
       ReplaceElemInGroups( elem, newElem, aMesh );
     }
     aMesh->RemoveElement( elem );
@@ -9225,6 +8995,7 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
                                              SMESH_MesherHelper& theHelper,
                                              const bool          theForce3d)
 {
+  //MESSAGE("convertElemToQuadratic");
   int nbElem = 0;
   if( !theSm ) return nbElem;
 
@@ -9250,18 +9021,20 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
       case SMDSEntity_Quad_Triangle:
       case SMDSEntity_Quad_Quadrangle:
       case SMDSEntity_Quad_Hexa:
+      case SMDSEntity_Quad_Penta:
         alreadyOK = !theHelper.GetIsBiQuadratic(); break;
 
       case SMDSEntity_BiQuad_Triangle:
       case SMDSEntity_BiQuad_Quadrangle:
       case SMDSEntity_TriQuad_Hexa:
+      case SMDSEntity_BiQuad_Penta:
         alreadyOK = theHelper.GetIsBiQuadratic();
         hasCentralNodes = true;
         break;
       default:
         alreadyOK = true;
       }
-      // take into account already present modium nodes
+      // take into account already present medium nodes
       switch ( aType ) {
       case SMDSAbs_Volume:
         theHelper.AddTLinks( static_cast< const SMDS_MeshVolume* >( elem )); break;
@@ -9280,14 +9053,14 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
     const int nbNodes = elem->NbCornerNodes();
     nodes.assign(elem->begin_nodes(), elem->end_nodes());
     if ( aGeomType == SMDSEntity_Polyhedra )
-      nbNodeInFaces = static_cast<const SMDS_VtkVolume* >( elem )->GetQuantities();
+      nbNodeInFaces = static_cast<const SMDS_MeshVolume* >( elem )->GetQuantities();
     else if ( aGeomType == SMDSEntity_Hexagonal_Prism )
       volumeToPolyhedron( elem, nodes, nbNodeInFaces );
 
     // remove a linear element
     GetMeshDS()->RemoveFreeElement(elem, theSm, /*fromGroups=*/false);
 
-    // remove central nodes of biquadratic elements (biquad->quad convertion)
+    // remove central nodes of biquadratic elements (biquad->quad conversion)
     if ( hasCentralNodes )
       for ( size_t i = nbNodes * 2; i < nodes.size(); ++i )
         if ( nodes[i]->NbInverseElements() == 0 )
@@ -9298,50 +9071,52 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
     switch( aType )
     {
     case SMDSAbs_Edge :
-      {
-        NewElem = theHelper.AddEdge(nodes[0], nodes[1], id, theForce3d);
-        break;
-      }
+    {
+      NewElem = theHelper.AddEdge(nodes[0], nodes[1], id, theForce3d);
+      break;
+    }
     case SMDSAbs_Face :
+    {
+      switch(nbNodes)
       {
-        switch(nbNodes)
-        {
-        case 3:
-          NewElem = theHelper.AddFace(nodes[0], nodes[1], nodes[2], id, theForce3d);
-          break;
-        case 4:
-          NewElem = theHelper.AddFace(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
-          break;
-        default:
-          NewElem = theHelper.AddPolygonalFace(nodes, id, theForce3d);
-        }
+      case 3:
+        NewElem = theHelper.AddFace(nodes[0], nodes[1], nodes[2], id, theForce3d);
+        break;
+      case 4:
+        NewElem = theHelper.AddFace(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
         break;
+      default:
+        NewElem = theHelper.AddPolygonalFace(nodes, id, theForce3d);
       }
+      break;
+    }
     case SMDSAbs_Volume :
+    {
+      switch( aGeomType )
       {
-        switch( aGeomType )
-        {
-        case SMDSEntity_Tetra:
-          NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
-          break;
-        case SMDSEntity_Pyramid:
-          NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], id, theForce3d);
-          break;
-        case SMDSEntity_Penta:
-          NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], nodes[5], id, theForce3d);
-          break;
-        case SMDSEntity_Hexa:
-        case SMDSEntity_Quad_Hexa:
-        case SMDSEntity_TriQuad_Hexa:
-          NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
-                                        nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
-          break;
-        case SMDSEntity_Hexagonal_Prism:
-        default:
-          NewElem = theHelper.AddPolyhedralVolume(nodes, nbNodeInFaces, id, theForce3d);
-        }
+      case SMDSEntity_Tetra:
+        NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+        break;
+      case SMDSEntity_Pyramid:
+        NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], id, theForce3d);
+        break;
+      case SMDSEntity_Penta:
+      case SMDSEntity_Quad_Penta:
+      case SMDSEntity_BiQuad_Penta:
+        NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], nodes[5], id, theForce3d);
+        break;
+      case SMDSEntity_Hexa:
+      case SMDSEntity_Quad_Hexa:
+      case SMDSEntity_TriQuad_Hexa:
+        NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+                                      nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
         break;
+      case SMDSEntity_Hexagonal_Prism:
+      default:
+        NewElem = theHelper.AddPolyhedralVolume(nodes, nbNodeInFaces, id, theForce3d);
       }
+      break;
+    }
     default :
       continue;
     }
@@ -9358,6 +9133,7 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
 
 void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theToBiQuad)
 {
+  //MESSAGE("ConvertToQuadratic "<< theForce3d << " " << theToBiQuad);
   SMESHDS_Mesh* meshDS = GetMeshDS();
 
   SMESH_MesherHelper aHelper(*myMesh);
@@ -9485,6 +9261,8 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT
         {
         case SMDSEntity_Quad_Hexa:    alreadyOK = !theToBiQuad; break;
         case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; break;
+        case SMDSEntity_Quad_Penta:   alreadyOK = !theToBiQuad; break;
+        case SMDSEntity_BiQuad_Penta: alreadyOK = theToBiQuad; break;
         default:                      alreadyOK = true;
         }
         if ( alreadyOK )
@@ -9496,7 +9274,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT
       const int id = volume->GetID();
       vector<const SMDS_MeshNode *> nodes (volume->begin_nodes(), volume->end_nodes());
       if ( type == SMDSEntity_Polyhedra )
-        nbNodeInFaces = static_cast<const SMDS_VtkVolume* >(volume)->GetQuantities();
+        nbNodeInFaces = static_cast<const SMDS_MeshVolume* >(volume)->GetQuantities();
       else if ( type == SMDSEntity_Hexagonal_Prism )
         volumeToPolyhedron( volume, nodes, nbNodeInFaces );
 
@@ -9522,8 +9300,13 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT
                                       nodes[3], nodes[4], id, theForce3d);
         break;
       case SMDSEntity_Penta:
+      case SMDSEntity_Quad_Penta:
+      case SMDSEntity_BiQuad_Penta:
         NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
                                       nodes[3], nodes[4], nodes[5], id, theForce3d);
+        for ( size_t i = 15; i < nodes.size(); ++i ) // rm central nodes
+          if ( nodes[i]->NbInverseElements() == 0 )
+            GetMeshDS()->RemoveFreeNode( nodes[i], /*sm=*/0, /*fromGroups=*/true );
         break;
       case SMDSEntity_Hexagonal_Prism:
       default:
@@ -9711,7 +9494,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool        theForce3d,
     if( newElem && smDS )
       smDS->AddElement( newElem );
 
-     // remove central nodes
+    // remove central nodes
     for ( size_t i = nodes.size() - nbCentralNodes; i < nodes.size(); ++i )
       if ( nodes[i]->NbInverseElements() == 0 )
         meshDS->RemoveFreeNode( nodes[i], smDS, /*fromGroups=*/true );
@@ -9845,7 +9628,7 @@ void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements)
   }
 
   // replace given elements by linear ones
-  SMDS_ElemIteratorPtr elemIt = elemSetIterator( theElements );
+  SMDS_ElemIteratorPtr elemIt = SMESHUtils::elemSetIterator( theElements );
   removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 );
 
   // we need to convert remaining elements whose all medium nodes are in mediumNodeIDs
@@ -9897,7 +9680,7 @@ void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements)
       }
     }
   }
-  elemIt = elemSetIterator( moreElemsToConvert );
+  elemIt = SMESHUtils::elemSetIterator( moreElemsToConvert );
   removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 );
 }
 
@@ -9914,10 +9697,8 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
                                    const SMDS_MeshNode* theSecondNode1,
                                    const SMDS_MeshNode* theSecondNode2)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
-  MESSAGE ("::::SewSideElements()");
   if ( theSide1.size() != theSide2.size() )
     return SEW_DIFF_NB_OF_ELEMENTS;
 
@@ -9936,7 +9717,7 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
   //    face does not exist
 
   SMESHDS_Mesh* aMesh = GetMeshDS();
-  // TODO algoritm not OK with vtkUnstructuredGrid: 2 meshes can't share nodes
+  // TODO algorithm not OK with vtkUnstructuredGrid: 2 meshes can't share nodes
   //SMDS_Mesh aTmpFacesMesh; // try to use the same mesh
   TIDSortedElemSet             faceSet1, faceSet2;
   set<const SMDS_MeshElement*> volSet1,  volSet2;
@@ -10192,12 +9973,12 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
 
   if ( faceSet1.size() != faceSet2.size() ) {
     // delete temporary faces: they are in reverseElements of actual nodes
-//    SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
-//    while ( tmpFaceIt->more() )
-//      aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
-//    list<const SMDS_MeshElement* >::iterator tmpFaceIt = tempFaceList.begin();
-//    for (; tmpFaceIt !=tempFaceList.end(); ++tmpFaceIt)
-//      aMesh->RemoveElement(*tmpFaceIt);
+    //    SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
+    //    while ( tmpFaceIt->more() )
+    //      aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
+    //    list<const SMDS_MeshElement* >::iterator tmpFaceIt = tempFaceList.begin();
+    //    for (; tmpFaceIt !=tempFaceList.end(); ++tmpFaceIt)
+    //      aMesh->RemoveElement(*tmpFaceIt);
     MESSAGE("Diff nb of faces");
     return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
   }
@@ -10255,20 +10036,8 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
         //cout << " F " << face[ iSide]->GetID() <<endl;
         faceSetPtr[ iSide ]->erase( face[ iSide ]);
         // put face nodes to fnodes
-        if ( face[ iSide ]->IsQuadratic() )
-        {
-          // use interlaced nodes iterator
-          const SMDS_VtkFace* F = dynamic_cast<const SMDS_VtkFace*>( face[ iSide ]);
-          if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
-          SMDS_ElemIteratorPtr nIter = F->interlacedNodesElemIterator();
-          while ( nIter->more() )
-            fnodes[ iSide ].push_back( cast2Node( nIter->next() ));
-        }
-        else
-        {
-          fnodes[ iSide ].assign( face[ iSide ]->begin_nodes(),
-                                  face[ iSide ]->end_nodes() );
-        }
+        SMDS_MeshElement::iterator nIt( face[ iSide ]->interlacedNodesIterator() ), nEnd;
+        fnodes[ iSide ].assign( nIt, nEnd );
         fnodes[ iSide ].push_back( fnodes[ iSide ].front());
       }
     }
@@ -10336,7 +10105,7 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
 
   if ( aResult == SEW_OK &&
        ( //linkIt[0] != linkList[0].end() ||
-         !faceSetPtr[0]->empty() || !faceSetPtr[1]->empty() )) {
+        !faceSetPtr[0]->empty() || !faceSetPtr[1]->empty() )) {
     MESSAGE( (linkIt[0] != linkList[0].end()) <<" "<< (faceSetPtr[0]->empty()) <<
              " " << (faceSetPtr[1]->empty()));
     aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
@@ -10347,9 +10116,9 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
   // ====================================================================
 
   // delete temporary faces
-//  SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
-//  while ( tmpFaceIt->more() )
-//    aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
+  //  SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
+  //  while ( tmpFaceIt->more() )
+  //    aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
   list<const SMDS_MeshElement* >::iterator tmpFaceIt = tempFaceList.begin();
   for (; tmpFaceIt !=tempFaceList.end(); ++tmpFaceIt)
     aMesh->RemoveElement(*tmpFaceIt);
@@ -10581,6 +10350,559 @@ SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
   return SEW_OK;
 }
 
+namespace // automatically find theAffectedElems for DoubleNodes()
+{
+  bool isOut( const SMDS_MeshNode* n, const gp_XYZ& norm, const SMDS_MeshElement* elem );
+
+  //--------------------------------------------------------------------------------
+  // Nodes shared by adjacent FissureBorder's.
+  // 1 node  if FissureBorder separates faces
+  // 2 nodes if FissureBorder separates volumes
+  struct SubBorder
+  {
+    const SMDS_MeshNode* _nodes[2];
+    int                  _nbNodes;
+
+    SubBorder( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2 = 0 )
+    {
+      _nodes[0] = n1;
+      _nodes[1] = n2;
+      _nbNodes = bool( n1 ) + bool( n2 );
+      if ( _nbNodes == 2 && n1 > n2 )
+        std::swap( _nodes[0], _nodes[1] );
+    }
+    bool operator<( const SubBorder& other ) const
+    {
+      for ( int i = 0; i < _nbNodes; ++i )
+      {
+        if ( _nodes[i] < other._nodes[i] ) return true;
+        if ( _nodes[i] > other._nodes[i] ) return false;
+      }
+      return false;
+    }
+  };
+
+  //--------------------------------------------------------------------------------
+  // Map a SubBorder to all FissureBorder it bounds
+  struct FissureBorder;
+  typedef std::map< SubBorder, std::vector< FissureBorder* > > TBorderLinks;
+  typedef TBorderLinks::iterator                               TMappedSub;
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Element border (volume facet or face edge) at a fissure
+   */
+  struct FissureBorder
+  {
+    std::vector< const SMDS_MeshNode* > _nodes;    // border nodes
+    const SMDS_MeshElement*             _elems[2]; // volume or face adjacent to fissure
+
+    std::vector< TMappedSub           > _mappedSubs;  // Sub() in TBorderLinks map
+    std::vector< const SMDS_MeshNode* > _sortedNodes; // to compare FissureBorder's
+
+    FissureBorder( FissureBorder && from ) // move constructor
+    {
+      std::swap( _nodes,       from._nodes );
+      std::swap( _sortedNodes, from._sortedNodes );
+      _elems[0] = from._elems[0];
+      _elems[1] = from._elems[1];
+    }
+
+    FissureBorder( const SMDS_MeshElement*                  elemToDuplicate,
+                   std::vector< const SMDS_MeshElement* > & adjElems)
+      : _nodes( elemToDuplicate->NbCornerNodes() )
+    {
+      for ( size_t i = 0; i < _nodes.size(); ++i )
+        _nodes[i] = elemToDuplicate->GetNode( i );
+
+      SMDSAbs_ElementType type = SMDSAbs_ElementType( elemToDuplicate->GetType() + 1 );
+      findAdjacent( type, adjElems );
+    }
+
+    FissureBorder( const SMDS_MeshNode**                    nodes,
+                   const size_t                             nbNodes,
+                   const SMDSAbs_ElementType                adjElemsType,
+                   std::vector< const SMDS_MeshElement* > & adjElems)
+      : _nodes( nodes, nodes + nbNodes )
+    {
+      findAdjacent( adjElemsType, adjElems );
+    }
+
+    void findAdjacent( const SMDSAbs_ElementType                adjElemsType,
+                       std::vector< const SMDS_MeshElement* > & adjElems)
+    {
+      _elems[0] = _elems[1] = 0;
+      adjElems.clear();
+      if ( SMDS_Mesh::GetElementsByNodes( _nodes, adjElems, adjElemsType ))
+        for ( size_t i = 0; i < adjElems.size() && i < 2; ++i )
+          _elems[i] = adjElems[i];
+    }
+
+    bool operator<( const FissureBorder& other ) const
+    {
+      return GetSortedNodes() < other.GetSortedNodes();
+    }
+
+    const std::vector< const SMDS_MeshNode* >& GetSortedNodes() const
+    {
+      if ( _sortedNodes.empty() && !_nodes.empty() )
+      {
+        FissureBorder* me = const_cast<FissureBorder*>( this );
+        me->_sortedNodes = me->_nodes;
+        std::sort( me->_sortedNodes.begin(), me->_sortedNodes.end() );
+      }
+      return _sortedNodes;
+    }
+
+    size_t NbSub() const
+    {
+      return _nodes.size();
+    }
+
+    SubBorder Sub(size_t i) const
+    {
+      return SubBorder( _nodes[i], NbSub() > 2 ? _nodes[ (i+1)%NbSub() ] : 0 );
+    }
+
+    void AddSelfTo( TBorderLinks& borderLinks )
+    {
+      _mappedSubs.resize( NbSub() );
+      for ( size_t i = 0; i < NbSub(); ++i )
+      {
+        TBorderLinks::iterator s2b =
+          borderLinks.insert( std::make_pair( Sub(i), TBorderLinks::mapped_type() )).first;
+        s2b->second.push_back( this );
+        _mappedSubs[ i ] = s2b;
+      }
+    }
+
+    void Clear()
+    {
+      _nodes.clear();
+    }
+
+    const SMDS_MeshElement* GetMarkedElem() const
+    {
+      if ( _nodes.empty() ) return 0; // cleared
+      if ( _elems[0] && _elems[0]->isMarked() ) return _elems[0];
+      if ( _elems[1] && _elems[1]->isMarked() ) return _elems[1];
+      return 0;
+    }
+
+    gp_XYZ GetNorm() const // normal to the border
+    {
+      gp_XYZ norm;
+      if ( _nodes.size() == 2 )
+      {
+        gp_XYZ avgNorm( 0,0,0 ); // sum of normals of adjacent faces
+        if ( SMESH_MeshAlgos::FaceNormal( _elems[0], norm ))
+          avgNorm += norm;
+        if ( SMESH_MeshAlgos::FaceNormal( _elems[1], norm ))
+          avgNorm += norm;
+
+        gp_XYZ bordDir( SMESH_NodeXYZ( _nodes[0] ) - SMESH_NodeXYZ( _nodes[1] ));
+        norm = bordDir ^ avgNorm;
+      }
+      else
+      {
+        SMESH_NodeXYZ p0( _nodes[0] );
+        SMESH_NodeXYZ p1( _nodes[1] );
+        SMESH_NodeXYZ p2( _nodes[2] );
+        norm = ( p0 - p1 ) ^ ( p2 - p1 );
+      }
+      if ( isOut( _nodes[0], norm, GetMarkedElem() ))
+        norm.Reverse();
+
+      return norm;
+    }
+
+    void ChooseSide() // mark an _elem located at positive side of fissure
+    {
+      _elems[0]->setIsMarked( true );
+      gp_XYZ norm = GetNorm();
+      double maxX = norm.Coord(1);
+      if ( Abs( maxX ) < Abs( norm.Coord(2)) ) maxX = norm.Coord(2);
+      if ( Abs( maxX ) < Abs( norm.Coord(3)) ) maxX = norm.Coord(3);
+      if ( maxX < 0 )
+      {
+        _elems[0]->setIsMarked( false );
+        _elems[1]->setIsMarked( true );
+      }
+    }
+
+  }; // struct FissureBorder
+
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Classifier of elements at fissure edge
+   */
+  class FissureNormal
+  {
+    std::vector< gp_XYZ > _normals;
+    bool                  _bothIn;
+
+  public:
+    void Add( const SMDS_MeshNode* n, const FissureBorder& bord )
+    {
+      _bothIn = false;
+      _normals.reserve(2);
+      _normals.push_back( bord.GetNorm() );
+      if ( _normals.size() == 2 )
+        _bothIn = !isOut( n, _normals[0], bord.GetMarkedElem() );
+    }
+
+    bool IsIn( const SMDS_MeshNode* n, const SMDS_MeshElement* elem ) const
+    {
+      bool isIn = false;
+      switch ( _normals.size() ) {
+      case 1:
+      {
+        isIn = !isOut( n, _normals[0], elem );
+        break;
+      }
+      case 2:
+      {
+        bool in1 = !isOut( n, _normals[0], elem );
+        bool in2 = !isOut( n, _normals[1], elem );
+        isIn = _bothIn ? ( in1 && in2 ) : ( in1 || in2 );
+      }
+      }
+      return isIn;
+    }
+  };
+
+  //================================================================================
+  /*!
+   * \brief Classify an element by a plane passing through a node
+   */
+  //================================================================================
+
+  bool isOut( const SMDS_MeshNode* n, const gp_XYZ& norm, const SMDS_MeshElement* elem )
+  {
+    SMESH_NodeXYZ p = n;
+    double sumDot = 0;
+    for ( int i = 0, nb = elem->NbCornerNodes(); i < nb; ++i )
+    {
+      SMESH_NodeXYZ pi = elem->GetNode( i );
+      sumDot += norm * ( pi - p );
+    }
+    return sumDot < -1e-100;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Find FissureBorder's by nodes to duplicate
+   */
+  //================================================================================
+
+  void findFissureBorders( const TIDSortedElemSet&        theNodes,
+                           std::vector< FissureBorder > & theFissureBorders )
+  {
+    TIDSortedElemSet::const_iterator nIt = theNodes.begin();
+    const SMDS_MeshNode* n = dynamic_cast< const SMDS_MeshNode*>( *nIt );
+    if ( !n ) return;
+    SMDSAbs_ElementType elemType = SMDSAbs_Volume;
+    if ( n->NbInverseElements( elemType ) == 0 )
+    {
+      elemType = SMDSAbs_Face;
+      if ( n->NbInverseElements( elemType ) == 0 )
+        return;
+    }
+    // unmark elements touching the fissure
+    for ( ; nIt != theNodes.end(); ++nIt )
+      SMESH_MeshAlgos::MarkElems( cast2Node(*nIt)->GetInverseElementIterator(), false );
+
+    // loop on elements touching the fissure to get their borders belonging to the fissure
+    std::set< FissureBorder >              fissureBorders;
+    std::vector< const SMDS_MeshElement* > adjElems;
+    std::vector< const SMDS_MeshNode* >    nodes;
+    SMDS_VolumeTool volTool;
+    for ( nIt = theNodes.begin(); nIt != theNodes.end(); ++nIt )
+    {
+      SMDS_ElemIteratorPtr invIt = cast2Node(*nIt)->GetInverseElementIterator( elemType );
+      while ( invIt->more() )
+      {
+        const SMDS_MeshElement* eInv = invIt->next();
+        if ( eInv->isMarked() ) continue;
+        eInv->setIsMarked( true );
+
+        if ( elemType == SMDSAbs_Volume )
+        {
+          volTool.Set( eInv );
+          int iQuad = eInv->IsQuadratic() ? 2 : 1;
+          for ( int iF = 0, nbF = volTool.NbFaces(); iF < nbF; ++iF )
+          {
+            const SMDS_MeshNode** nn = volTool.GetFaceNodes( iF );
+            int                  nbN = volTool.NbFaceNodes( iF ) / iQuad;
+            nodes.clear();
+            bool allOnFissure = true;
+            for ( int iN = 0; iN < nbN  && allOnFissure; iN += iQuad )
+              if (( allOnFissure = theNodes.count( nn[ iN ])))
+                nodes.push_back( nn[ iN ]);
+            if ( allOnFissure )
+              fissureBorders.insert( std::move( FissureBorder( &nodes[0], nodes.size(),
+                                                               elemType, adjElems )));
+          }
+        }
+        else // elemType == SMDSAbs_Face
+        {
+          const SMDS_MeshNode* nn[2] = { eInv->GetNode( eInv->NbCornerNodes()-1 ), 0 };
+          bool            onFissure0 = theNodes.count( nn[0] ), onFissure1;
+          for ( int iN = 0, nbN = eInv->NbCornerNodes(); iN < nbN; ++iN )
+          {
+            nn[1]      = eInv->GetNode( iN );
+            onFissure1 = theNodes.count( nn[1] );
+            if ( onFissure0 && onFissure1 )
+              fissureBorders.insert( std::move( FissureBorder( nn, 2, elemType, adjElems )));
+            nn[0]      = nn[1];
+            onFissure0 = onFissure1;
+          }
+        }
+      }
+    }
+
+    theFissureBorders.reserve( theFissureBorders.size() + fissureBorders.size());
+    std::set< FissureBorder >::iterator bord = fissureBorders.begin();
+    for ( ; bord != fissureBorders.end(); ++bord )
+    {
+      theFissureBorders.push_back( std::move( const_cast<FissureBorder&>( *bord ) ));
+    }
+    return;
+  } // findFissureBorders()
+
+  //================================================================================
+  /*!
+   * \brief Find elements on one side of a fissure defined by elements or nodes to duplicate
+   *  \param [in] theElemsOrNodes - elements or nodes to duplicate
+   *  \param [in] theNodesNot - nodes not to duplicate
+   *  \param [out] theAffectedElems - the found elements
+   */
+  //================================================================================
+
+  void findAffectedElems( const TIDSortedElemSet& theElemsOrNodes,
+                          TIDSortedElemSet&       theAffectedElems)
+  {
+    if ( theElemsOrNodes.empty() ) return;
+
+    // find FissureBorder's
+
+    std::vector< FissureBorder >           fissure;
+    std::vector< const SMDS_MeshElement* > elemsByFacet;
+
+    TIDSortedElemSet::const_iterator elIt = theElemsOrNodes.begin();
+    if ( (*elIt)->GetType() == SMDSAbs_Node )
+    {
+      findFissureBorders( theElemsOrNodes, fissure );
+    }
+    else
+    {
+      fissure.reserve( theElemsOrNodes.size() );
+      for ( ; elIt != theElemsOrNodes.end(); ++elIt )
+        fissure.push_back( std::move( FissureBorder( *elIt, elemsByFacet )));
+    }
+    if ( fissure.empty() )
+      return;
+
+    // fill borderLinks
+
+    TBorderLinks borderLinks;
+
+    for ( size_t i = 0; i < fissure.size(); ++i )
+    {
+      fissure[i].AddSelfTo( borderLinks );
+    }
+
+    // get theAffectedElems
+
+    // unmark elements having nodes on the fissure, theAffectedElems elements will be marked
+    for ( size_t i = 0; i < fissure.size(); ++i )
+      for ( size_t j = 0; j < fissure[i]._nodes.size(); ++j )
+      {
+        SMESH_MeshAlgos::MarkElemNodes( fissure[i]._nodes[j]->GetInverseElementIterator(),
+                                        false, /*markElem=*/true );
+      }
+
+    std::vector<const SMDS_MeshNode *>                 facetNodes;
+    std::map< const SMDS_MeshNode*, FissureNormal >    fissEdgeNodes2Norm;
+    boost::container::flat_set< const SMDS_MeshNode* > fissureNodes;
+
+    // choose a side of fissure
+    fissure[0].ChooseSide();
+    theAffectedElems.insert( fissure[0].GetMarkedElem() );
+
+    size_t nbCheckedBorders = 0;
+    while ( nbCheckedBorders < fissure.size() )
+    {
+      // find a FissureBorder to treat
+      FissureBorder* bord = 0;
+      for ( size_t i = 0; i < fissure.size()  && !bord; ++i )
+        if ( fissure[i].GetMarkedElem() )
+          bord = & fissure[i];
+      for ( size_t i = 0; i < fissure.size()  && !bord; ++i )
+        if ( fissure[i].NbSub() > 0 && fissure[i]._elems[0] )
+        {
+          bord = & fissure[i];
+          bord->ChooseSide();
+          theAffectedElems.insert( bord->GetMarkedElem() );
+        }
+      if ( !bord ) return;
+      ++nbCheckedBorders;
+
+      // treat FissureBorder's linked to bord
+      fissureNodes.clear();
+      fissureNodes.insert( bord->_nodes.begin(), bord->_nodes.end() );
+      for ( size_t i = 0; i < bord->NbSub(); ++i )
+      {
+        TBorderLinks::iterator l2b = bord->_mappedSubs[ i ];
+        if ( l2b == borderLinks.end() || l2b->second.empty() ) continue;
+        std::vector< FissureBorder* >& linkedBorders = l2b->second;
+        const SubBorder&                          sb = l2b->first;
+        const SMDS_MeshElement*             bordElem = bord->GetMarkedElem();
+
+        if ( linkedBorders.size() == 1 ) // fissure edge reached, fill fissEdgeNodes2Norm
+        {
+          for ( int j = 0; j < sb._nbNodes; ++j )
+            fissEdgeNodes2Norm[ sb._nodes[j] ].Add( sb._nodes[j], *bord );
+          continue;
+        }
+
+        // add to theAffectedElems elems sharing nodes of a SubBorder and a node of bordElem
+        // until an elem adjacent to a neighbour FissureBorder is found
+        facetNodes.clear();
+        facetNodes.insert( facetNodes.end(), sb._nodes, sb._nodes + sb._nbNodes );
+        facetNodes.resize( sb._nbNodes + 1 );
+
+        while ( bordElem )
+        {
+          // check if bordElem is adjacent to a neighbour FissureBorder
+          for ( size_t j = 0; j < linkedBorders.size(); ++j )
+          {
+            FissureBorder* bord2 = linkedBorders[j];
+            if ( bord2 == bord ) continue;
+            if ( bordElem == bord2->_elems[0] || bordElem == bord2->_elems[1] )
+              bordElem = 0;
+            else
+              fissureNodes.insert( bord2->_nodes.begin(), bord2->_nodes.end() );
+          }
+          if ( !bordElem )
+            break;
+
+          // find the next bordElem
+          const SMDS_MeshElement* nextBordElem = 0;
+          for ( int iN = 0, nbN = bordElem->NbCornerNodes(); iN < nbN  && !nextBordElem; ++iN )
+          {
+            const SMDS_MeshNode* n = bordElem->GetNode( iN );
+            if ( fissureNodes.count( n )) continue;
+
+            facetNodes[ sb._nbNodes ] = n;
+            elemsByFacet.clear();
+            if ( SMDS_Mesh::GetElementsByNodes( facetNodes, elemsByFacet ) > 1 )
+            {
+              for ( size_t iE = 0; iE < elemsByFacet.size(); ++iE )
+                if ( elemsByFacet[ iE ] != bordElem &&
+                     !elemsByFacet[ iE ]->isMarked() )
+                {
+                  theAffectedElems.insert( elemsByFacet[ iE ]);
+                  elemsByFacet[ iE ]->setIsMarked( true );
+                  if ( elemsByFacet[ iE ]->GetType() == bordElem->GetType() )
+                    nextBordElem = elemsByFacet[ iE ];
+                }
+            }
+          }
+          bordElem = nextBordElem;
+
+        } // while ( bordElem )
+
+        linkedBorders.clear(); // not to treat this link any more
+
+      } // loop on SubBorder's of a FissureBorder
+
+      bord->Clear();
+
+    } // loop on FissureBorder's
+
+
+    // add elements sharing only one node of the fissure, except those sharing fissure edge nodes
+
+    // mark nodes of theAffectedElems
+    SMESH_MeshAlgos::MarkElemNodes( theAffectedElems.begin(), theAffectedElems.end(), true );
+
+    // unmark nodes of the fissure
+    elIt = theElemsOrNodes.begin();
+    if ( (*elIt)->GetType() == SMDSAbs_Node )
+      SMESH_MeshAlgos::MarkElems( elIt, theElemsOrNodes.end(), false );
+    else
+      SMESH_MeshAlgos::MarkElemNodes( elIt, theElemsOrNodes.end(), false );
+
+    std::vector< gp_XYZ > normVec;
+
+    // loop on nodes of the fissure, add elements having marked nodes
+    for ( elIt = theElemsOrNodes.begin(); elIt != theElemsOrNodes.end(); ++elIt )
+    {
+      const SMDS_MeshElement* e = (*elIt);
+      if ( e->GetType() != SMDSAbs_Node )
+        e->setIsMarked( true ); // avoid adding a fissure element
+
+      for ( int iN = 0, nbN = e->NbCornerNodes(); iN < nbN; ++iN )
+      {
+        const SMDS_MeshNode* n = e->GetNode( iN );
+        if ( fissEdgeNodes2Norm.count( n ))
+          continue;
+
+        SMDS_ElemIteratorPtr invIt = n->GetInverseElementIterator();
+        while ( invIt->more() )
+        {
+          const SMDS_MeshElement* eInv = invIt->next();
+          if ( eInv->isMarked() ) continue;
+          eInv->setIsMarked( true );
+
+          SMDS_ElemIteratorPtr nIt = eInv->nodesIterator();
+          while( nIt->more() )
+            if ( nIt->next()->isMarked())
+            {
+              theAffectedElems.insert( eInv );
+              SMESH_MeshAlgos::MarkElems( eInv->nodesIterator(), true );
+              n->setIsMarked( false );
+              break;
+            }
+        }
+      }
+    }
+
+    // add elements on the fissure edge
+    std::map< const SMDS_MeshNode*, FissureNormal >::iterator n2N;
+    for ( n2N = fissEdgeNodes2Norm.begin(); n2N != fissEdgeNodes2Norm.end(); ++n2N )
+    {
+      const SMDS_MeshNode* edgeNode = n2N->first;
+      const FissureNormal & normals = n2N->second;
+
+      SMDS_ElemIteratorPtr invIt = edgeNode->GetInverseElementIterator();
+      while ( invIt->more() )
+      {
+        const SMDS_MeshElement* eInv = invIt->next();
+        if ( eInv->isMarked() ) continue;
+        eInv->setIsMarked( true );
+
+        // classify eInv using normals
+        bool toAdd = normals.IsIn( edgeNode, eInv );
+        if ( toAdd ) // check if all nodes lie on the fissure edge
+        {
+          bool notOnEdge = false;
+          for ( int iN = 0, nbN = eInv->NbCornerNodes(); iN < nbN  && !notOnEdge; ++iN )
+            notOnEdge = !fissEdgeNodes2Norm.count( eInv->GetNode( iN ));
+          toAdd = notOnEdge;
+        }
+        if ( toAdd )
+        {
+          theAffectedElems.insert( eInv );
+        }
+      }
+    }
+
+    return;
+  } // findAffectedElems()
+} // namespace
+
 //================================================================================
 /*!
  * \brief Create elements equal (on same nodes) to given ones
@@ -10596,9 +10918,8 @@ void SMESH_MeshEditor::DoubleElements( const TIDSortedElemSet& theElements )
 
   // get an element type and an iterator over elements
 
-  SMDSAbs_ElementType type;
+  SMDSAbs_ElementType type = SMDSAbs_All;
   SMDS_ElemIteratorPtr elemIt;
-  vector< const SMDS_MeshElement* > allElems;
   if ( theElements.empty() )
   {
     if ( mesh->NbNodes() == 0 )
@@ -10612,21 +10933,19 @@ void SMESH_MeshEditor::DoubleElements( const TIDSortedElemSet& theElements )
       if ( mesh->GetMeshInfo().NbElements( types[i] ))
       {
         type = types[i];
+        elemIt = mesh->elementsIterator( type );
         break;
       }
-    // put all elements in the vector <allElems>
-    allElems.reserve( mesh->GetMeshInfo().NbElements( type ));
-    elemIt = mesh->elementsIterator( type );
-    while ( elemIt->more() )
-      allElems.push_back( elemIt->next());
-    elemIt = elemSetIterator( allElems );
   }
   else
   {
     type = (*theElements.begin())->GetType();
-    elemIt = elemSetIterator( theElements );
+    elemIt = SMESHUtils::elemSetIterator( theElements );
   }
 
+  // un-mark all elements to avoid duplicating just created elements
+  SMESH_MeshAlgos::MarkElems( mesh->elementsIterator( type ), false );
+
   // duplicate elements
 
   ElemFeatures elemType;
@@ -10635,13 +10954,14 @@ void SMESH_MeshEditor::DoubleElements( const TIDSortedElemSet& theElements )
   while ( elemIt->more() )
   {
     const SMDS_MeshElement* elem = elemIt->next();
-    if ( elem->GetType() != type )
+    if ( elem->GetType() != type || elem->isMarked() )
       continue;
 
     elemType.Init( elem, /*basicOnly=*/false );
     nodes.assign( elem->begin_nodes(), elem->end_nodes() );
 
-    AddElement( nodes, elemType );
+    if ( const SMDS_MeshElement* newElem = AddElement( nodes, elemType ))
+      newElem->setIsMarked( true );
   }
 }
 
@@ -10661,8 +10981,7 @@ bool SMESH_MeshEditor::DoubleNodes( const TIDSortedElemSet& theElems,
                                     const TIDSortedElemSet& theNodesNot,
                                     const TIDSortedElemSet& theAffectedElems )
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
   if ( theElems.size() == 0 )
     return false;
@@ -10698,7 +11017,6 @@ bool SMESH_MeshEditor::doubleNodes(SMESHDS_Mesh*           theMeshDS,
                                    TNodeNodeMap&           theNodeNodeMap,
                                    const bool              theIsDoubleElem )
 {
-  MESSAGE("doubleNodes");
   // iterate through element and duplicate them (by nodes duplication)
   bool res = false;
   std::vector<const SMDS_MeshNode*> newNodes;
@@ -10708,8 +11026,8 @@ bool SMESH_MeshEditor::doubleNodes(SMESHDS_Mesh*           theMeshDS,
   for ( ;  elemItr != theElems.end(); ++elemItr )
   {
     const SMDS_MeshElement* anElem = *elemItr;
-    if (!anElem)
-      continue;
+    // if (!anElem)
+    //   continue;
 
     // duplicate nodes to duplicate element
     bool isDuplicate = false;
@@ -10731,7 +11049,7 @@ bool SMESH_MeshEditor::doubleNodes(SMESHDS_Mesh*           theMeshDS,
         aNewNode = theMeshDS->AddNode( aCurrNode->X(), aCurrNode->Y(), aCurrNode->Z() );
         copyPosition( aCurrNode, aNewNode );
         theNodeNodeMap[ aCurrNode ] = aNewNode;
-        myLastCreatedNodes.Append( aNewNode );
+        myLastCreatedNodes.push_back( aNewNode );
       }
       isDuplicate |= (aCurrNode != aNewNode);
       newNodes[ ind++ ] = aNewNode;
@@ -10763,9 +11081,7 @@ bool SMESH_MeshEditor::doubleNodes(SMESHDS_Mesh*           theMeshDS,
 bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes,
                                     const std::list< int >& theListOfModifiedElems )
 {
-  MESSAGE("DoubleNodes");
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
+  ClearLastCreated();
 
   if ( theListOfNodes.size() == 0 )
     return false;
@@ -10781,8 +11097,7 @@ bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes,
   std::list< int >::const_iterator aNodeIter;
   for ( aNodeIter = theListOfNodes.begin(); aNodeIter != theListOfNodes.end(); ++aNodeIter )
   {
-    int aCurr = *aNodeIter;
-    SMDS_MeshNode* aNode = (SMDS_MeshNode*)aMeshDS->FindNode( aCurr );
+    const SMDS_MeshNode* aNode = aMeshDS->FindNode( *aNodeIter );
     if ( !aNode )
       continue;
 
@@ -10793,54 +11108,32 @@ bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes,
     {
       copyPosition( aNode, aNewNode );
       anOldNodeToNewNode[ aNode ] = aNewNode;
-      myLastCreatedNodes.Append( aNewNode );
+      myLastCreatedNodes.push_back( aNewNode );
     }
   }
 
-  // Create map of new nodes for modified elements
+  // Change nodes of elements
 
-  std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> > anElemToNodes;
+  std::vector<const SMDS_MeshNode*> aNodeArr;
 
   std::list< int >::const_iterator anElemIter;
-  for ( anElemIter = theListOfModifiedElems.begin();
-        anElemIter != theListOfModifiedElems.end(); ++anElemIter )
+  for ( anElemIter =  theListOfModifiedElems.begin();
+        anElemIter != theListOfModifiedElems.end();
+        anElemIter++ )
   {
-    int aCurr = *anElemIter;
-    SMDS_MeshElement* anElem = (SMDS_MeshElement*)aMeshDS->FindElement( aCurr );
+    const SMDS_MeshElement* anElem = aMeshDS->FindElement( *anElemIter );
     if ( !anElem )
       continue;
 
-    vector<const SMDS_MeshNode*> aNodeArr( anElem->NbNodes() );
-
-    SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
-    int ind = 0;
-    while ( anIter->more() )
+    aNodeArr.assign( anElem->begin_nodes(), anElem->end_nodes() );
+    for( size_t i = 0; i < aNodeArr.size(); ++i )
     {
-      SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next();
-      if ( aCurr && anOldNodeToNewNode.find( aCurrNode ) != anOldNodeToNewNode.end() )
-      {
-        const SMDS_MeshNode* aNewNode = anOldNodeToNewNode[ aCurrNode ];
-        aNodeArr[ ind++ ] = aNewNode;
-      }
-      else
-        aNodeArr[ ind++ ] = aCurrNode;
+      std::map< const SMDS_MeshNode*, const SMDS_MeshNode* >::iterator n2n =
+        anOldNodeToNewNode.find( aNodeArr[ i ]);
+      if ( n2n != anOldNodeToNewNode.end() )
+        aNodeArr[ i ] = n2n->second;
     }
-    anElemToNodes[ anElem ] = aNodeArr;
-  }
-
-  // Change nodes of elements
-
-  std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> >::iterator
-    anElemToNodesIter = anElemToNodes.begin();
-  for ( ; anElemToNodesIter != anElemToNodes.end(); ++anElemToNodesIter )
-  {
-    const SMDS_MeshElement* anElem = anElemToNodesIter->first;
-    vector<const SMDS_MeshNode*> aNodeArr = anElemToNodesIter->second;
-    if ( anElem )
-      {
-      MESSAGE("ChangeElementNodes");
-      aMeshDS->ChangeElementNodes( anElem, &aNodeArr[ 0 ], anElem->NbNodes() );
-      }
+    aMeshDS->ChangeElementNodes( anElem, &aNodeArr[ 0 ], aNodeArr.size() );
   }
 
   return true;
@@ -10850,8 +11143,8 @@ namespace {
 
   //================================================================================
   /*!
-  \brief Check if element located inside shape
-  \return TRUE if IN or ON shape, FALSE otherwise
+    \brief Check if element located inside shape
+    \return TRUE if IN or ON shape, FALSE otherwise
   */
   //================================================================================
 
@@ -10861,9 +11154,8 @@ namespace {
                 const double            theTol)
   {
     gp_XYZ centerXYZ (0, 0, 0);
-    SMDS_ElemIteratorPtr aNodeItr = theElem->nodesIterator();
-    while (aNodeItr->more())
-      centerXYZ += SMESH_TNodeXYZ(cast2Node( aNodeItr->next()));
+    for ( SMDS_ElemIteratorPtr aNodeItr = theElem->nodesIterator(); aNodeItr->more(); )
+      centerXYZ += SMESH_NodeXYZ( aNodeItr->next() );
 
     gp_Pnt aPnt = centerXYZ / theElem->NbNodes();
     theClassifier.Perform(aPnt, theTol);
@@ -10918,9 +11210,9 @@ namespace {
          (select elements with a gravity center on the side given by faces normals).
          This mode (null shape) is faster, but works only when theElems are faces, with coherents orientations.
          The replicated nodes should be associated to affected elements.
-  \return groups of affected elements
+  \return true
   \sa DoubleNodeElemGroupsInRegion()
- */
+*/
 //================================================================================
 
 bool SMESH_MeshEditor::AffectedElemGroupsInRegion( const TIDSortedElemSet& theElems,
@@ -10930,101 +11222,7 @@ bool SMESH_MeshEditor::AffectedElemGroupsInRegion( const TIDSortedElemSet& theEl
 {
   if ( theShape.IsNull() )
   {
-    std::set<const SMDS_MeshNode*> alreadyCheckedNodes;
-    std::set<const SMDS_MeshElement*> alreadyCheckedElems;
-    std::set<const SMDS_MeshElement*> edgesToCheck;
-    alreadyCheckedNodes.clear();
-    alreadyCheckedElems.clear();
-    edgesToCheck.clear();
-
-    // --- iterates on elements to be replicated and get elements by back references from their nodes
-
-    TIDSortedElemSet::const_iterator elemItr = theElems.begin();
-    int ielem;
-    for ( ielem=1;  elemItr != theElems.end(); ++elemItr )
-    {
-      SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr;
-      if (!anElem || (anElem->GetType() != SMDSAbs_Face))
-        continue;
-      gp_XYZ normal;
-      SMESH_MeshAlgos::FaceNormal( anElem, normal, /*normalized=*/true );
-      MESSAGE("element " << ielem++ <<  " normal " << normal.X() << " " << normal.Y() << " " << normal.Z());
-      std::set<const SMDS_MeshNode*> nodesElem;
-      nodesElem.clear();
-      SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator();
-      while ( nodeItr->more() )
-      {
-        const SMDS_MeshNode* aNode = cast2Node(nodeItr->next());
-        nodesElem.insert(aNode);
-      }
-      std::set<const SMDS_MeshNode*>::iterator nodit = nodesElem.begin();
-      for (; nodit != nodesElem.end(); nodit++)
-      {
-        MESSAGE("  noeud ");
-        const SMDS_MeshNode* aNode = *nodit;
-        if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() )
-          continue;
-        if (alreadyCheckedNodes.find(aNode) != alreadyCheckedNodes.end())
-          continue;
-        alreadyCheckedNodes.insert(aNode);
-        SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator();
-        while ( backElemItr->more() )
-        {
-          MESSAGE("    backelem ");
-          const SMDS_MeshElement* curElem = backElemItr->next();
-          if (alreadyCheckedElems.find(curElem) != alreadyCheckedElems.end())
-            continue;
-          if (theElems.find(curElem) != theElems.end())
-            continue;
-          alreadyCheckedElems.insert(curElem);
-          double x=0, y=0, z=0;
-          int nb = 0;
-          SMDS_ElemIteratorPtr nodeItr2 = curElem->nodesIterator();
-          while ( nodeItr2->more() )
-          {
-            const SMDS_MeshNode* anotherNode = cast2Node(nodeItr2->next());
-            x += anotherNode->X();
-            y += anotherNode->Y();
-            z += anotherNode->Z();
-            nb++;
-          }
-          gp_XYZ p;
-          p.SetCoord( x/nb -aNode->X(),
-                      y/nb -aNode->Y(),
-                      z/nb -aNode->Z() );
-          MESSAGE("      check " << p.X() << " " << p.Y() << " " << p.Z());
-          if (normal*p > 0)
-          {
-            MESSAGE("    --- inserted")
-            theAffectedElems.insert( curElem );
-          }
-          else if (curElem->GetType() == SMDSAbs_Edge)
-            edgesToCheck.insert(curElem);
-        }
-      }
-    }
-    // --- add also edges lying on the set of faces (all nodes in alreadyCheckedNodes)
-    std::set<const SMDS_MeshElement*>::iterator eit = edgesToCheck.begin();
-    for( ; eit != edgesToCheck.end(); eit++)
-    {
-      bool onside = true;
-      const SMDS_MeshElement* anEdge = *eit;
-      SMDS_ElemIteratorPtr nodeItr = anEdge->nodesIterator();
-      while ( nodeItr->more() )
-      {
-        const SMDS_MeshNode* aNode = cast2Node(nodeItr->next());
-        if (alreadyCheckedNodes.find(aNode) == alreadyCheckedNodes.end())
-        {
-          onside = false;
-          break;
-        }
-      }
-      if (onside)
-      {
-        MESSAGE("    --- edge onside inserted")
-        theAffectedElems.insert(anEdge);
-      }
-    }
+    findAffectedElems( theElems, theAffectedElems );
   }
   else
   {
@@ -11043,29 +11241,23 @@ bool SMESH_MeshEditor::AffectedElemGroupsInRegion( const TIDSortedElemSet& theEl
 
     // iterates on indicated elements and get elements by back references from their nodes
     TIDSortedElemSet::const_iterator elemItr = theElems.begin();
-    int ielem;
-    for ( ielem = 1;  elemItr != theElems.end(); ++elemItr )
+    for ( ;  elemItr != theElems.end(); ++elemItr )
     {
-      MESSAGE("element " << ielem++);
-      SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr;
-      if (!anElem)
-        continue;
+      SMDS_MeshElement*     anElem = (SMDS_MeshElement*)*elemItr;
       SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator();
       while ( nodeItr->more() )
       {
-        MESSAGE("  noeud ");
         const SMDS_MeshNode* aNode = cast2Node(nodeItr->next());
         if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() )
           continue;
         SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator();
         while ( backElemItr->more() )
         {
-          MESSAGE("    backelem ");
           const SMDS_MeshElement* curElem = backElemItr->next();
           if ( curElem && theElems.find(curElem) == theElems.end() &&
-              ( bsc3d.get() ?
-                isInside( curElem, *bsc3d, aTol ) :
-                isInside( curElem, *aFaceClassifier, aTol )))
+               ( bsc3d.get() ?
+                 isInside( curElem, *bsc3d, aTol ) :
+                 isInside( curElem, *aFaceClassifier, aTol )))
             theAffectedElems.insert( curElem );
         }
       }
@@ -11146,10 +11338,6 @@ bool SMESH_MeshEditor::DoubleNodesInRegion( const TIDSortedElemSet& theElems,
  */
 double SMESH_MeshEditor::OrientedAngle(const gp_Pnt& p0, const gp_Pnt& p1, const gp_Pnt& g1, const gp_Pnt& g2)
 {
-//  MESSAGE("    p0: " << p0.X() << " " << p0.Y() << " " << p0.Z());
-//  MESSAGE("    p1: " << p1.X() << " " << p1.Y() << " " << p1.Z());
-//  MESSAGE("    g1: " << g1.X() << " " << g1.Y() << " " << g1.Z());
-//  MESSAGE("    g2: " << g2.X() << " " << g2.Y() << " " << g2.Z());
   gp_Vec vref(p0, p1);
   gp_Vec v1(p0, g1);
   gp_Vec v2(p0, g2);
@@ -11184,14 +11372,14 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
                                                      bool                                 createJointElems,
                                                      bool                                 onAllBoundaries)
 {
-  MESSAGE("----------------------------------------------");
-  MESSAGE("SMESH_MeshEditor::doubleNodesOnGroupBoundaries");
-  MESSAGE("----------------------------------------------");
+  // MESSAGE("----------------------------------------------");
+  // MESSAGE("SMESH_MeshEditor::doubleNodesOnGroupBoundaries");
+  // MESSAGE("----------------------------------------------");
 
   SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS();
   meshDS->BuildDownWardConnectivity(true);
   CHRONO(50);
-  SMDS_UnstructuredGrid *grid = meshDS->getGrid();
+  SMDS_UnstructuredGrid *grid = meshDS->GetGrid();
 
   // --- build the list of faces shared by 2 domains (group of elements), with their domain and volume indexes
   //     build the list of cells with only a node or an edge on the border, with their domain and volume indexes
@@ -11209,7 +11397,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   std::set<int> emptySet;
   emptyMap.clear();
 
-  MESSAGE(".. Number of domains :"<<theElems.size());
+  //MESSAGE(".. Number of domains :"<<theElems.size());
 
   TIDSortedElemSet theRestDomElems;
   const int iRestDom  = -1;
@@ -11249,7 +11437,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
     //     and corresponding volume of this domain, for each shared face.
     //     a volume has a face shared by 2 domains if it has a neighbor which is not in his domain.
 
-    MESSAGE("... Neighbors of domain #" << idom);
+    //MESSAGE("... Neighbors of domain #" << idom);
     const TIDSortedElemSet& domain = theElems[idom];
     TIDSortedElemSet::const_iterator elemItr = domain.begin();
     for (; elemItr != domain.end(); ++elemItr)
@@ -11257,7 +11445,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       const SMDS_MeshElement* anElem = *elemItr;
       if (!anElem)
         continue;
-      int vtkId = anElem->getVtkId();
+      int vtkId = anElem->GetVtkID();
       //MESSAGE("  vtkId " << vtkId << " smdsId " << anElem->GetID());
       int neighborsVtkIds[NBMAXNEIGHBORS];
       int downIds[NBMAXNEIGHBORS];
@@ -11265,7 +11453,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId);
       for (int n = 0; n < nbNeighbors; n++)
       {
-        int smdsId = meshDS->fromVtkToSmds(neighborsVtkIds[n]);
+        int smdsId = meshDS->FromVtkToSmds(neighborsVtkIds[n]);
         const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
         if (elem && ! domain.count(elem)) // neighbor is in another domain : face is shared
         {
@@ -11328,7 +11516,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
         for (int i=0; i<l.ncells; i++)
         {
           int vtkId = l.cells[i];
-          const SMDS_MeshElement* anElem = GetMeshDS()->FindElement(GetMeshDS()->fromVtkToSmds(vtkId));
+          const SMDS_MeshElement* anElem = GetMeshDS()->FindElement(GetMeshDS()->FromVtkToSmds(vtkId));
           if (!domain.count(anElem))
             continue;
           int vtkType = grid->GetCellType(vtkId);
@@ -11360,7 +11548,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   std::map<int, std::vector<int> > mutipleNodes; // nodes multi domains with domain order
   std::map<int, std::vector<int> > mutipleNodesToFace; // nodes multi domains with domain order to transform in Face (junction between 3 or more 2D domains)
 
-  MESSAGE(".. Duplication of the nodes");
+  //MESSAGE(".. Duplication of the nodes");
   for (int idomain = idom0; idomain < nbDomains; idomain++)
   {
     itface = faceDomains.begin();
@@ -11412,7 +11600,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
             double *coords = grid->GetPoint(oldId);
             SMDS_MeshNode *newNode = meshDS->AddNode(coords[0], coords[1], coords[2]);
             copyPosition( meshDS->FindNodeVtk( oldId ), newNode );
-            int newId = newNode->getVtkId();
+            int newId = newNode->GetVtkID();
             nodeDomains[oldId][idom] = newId; // cloned node for other domains
             //MESSAGE("-+-+-c     oldNode " << oldId << " domain " << idomain << " newNode " << newId << " domain " << idom << " size=" <<nodeDomains[oldId].size());
           }
@@ -11421,7 +11609,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
     }
   }
 
-  MESSAGE(".. Creation of elements");
+  //MESSAGE(".. Creation of elements");
   for (int idomain = idom0; idomain < nbDomains; idomain++)
   {
     itface = faceDomains.begin();
@@ -11485,7 +11673,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
                 gp_Pnt p1(coords[0], coords[1], coords[2]);
                 gp_Pnt gref;
                 int vtkVolIds[1000];  // an edge can belong to a lot of volumes
-                map<int, SMDS_VtkVolume*> domvol; // domain --> a volume with the edge
+                map<int, SMDS_MeshVolume*> domvol; // domain --> a volume with the edge
                 map<int, double> angleDom; // oriented angles between planes defined by edge and volume centers
                 int nbvol = grid->GetParentVolumes(vtkVolIds, downEdgeIds[ie], edgeType[ie]);
                 for ( size_t id = 0; id < doms.size(); id++ )
@@ -11494,26 +11682,31 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
                   const TIDSortedElemSet& domain = (idom == iRestDom) ? theRestDomElems : theElems[idom];
                   for ( int ivol = 0; ivol < nbvol; ivol++ )
                   {
-                    int smdsId = meshDS->fromVtkToSmds(vtkVolIds[ivol]);
-                    SMDS_MeshElement* elem = (SMDS_MeshElement*)meshDS->FindElement(smdsId);
+                    int smdsId = meshDS->FromVtkToSmds(vtkVolIds[ivol]);
+                    const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
                     if (domain.count(elem))
                     {
-                      SMDS_VtkVolume* svol = dynamic_cast<SMDS_VtkVolume*>(elem);
-                      domvol[idom] = svol;
+                      const SMDS_MeshVolume* svol = SMDS_Mesh::DownCast<SMDS_MeshVolume>(elem);
+                      domvol[idom] = (SMDS_MeshVolume*) svol;
                       //MESSAGE("  domain " << idom << " volume " << elem->GetID());
-                      double values[3];
+                      double values[3] = { 0,0,0 };
                       vtkIdType npts = 0;
                       vtkIdType* pts = 0;
                       grid->GetCellPoints(vtkVolIds[ivol], npts, pts);
-                      SMDS_VtkVolume::gravityCenter(grid, pts, npts, values);
-                      if (id ==0)
+                      for ( vtkIdType i = 0; i < npts; ++i )
+                      {
+                        double *coords = grid->GetPoint( pts[i] );
+                        for ( int j = 0; j < 3; ++j )
+                          values[j] += coords[j] / npts;
+                      }
+                      if ( id == 0 )
                       {
-                        gref.SetXYZ(gp_XYZ(values[0], values[1], values[2]));
+                        gref.SetCoord( values[0], values[1], values[2] );
                         angleDom[idom] = 0;
                       }
                       else
                       {
-                        gp_Pnt g(values[0], values[1], values[2]);
+                        gp_Pnt g( values[0], values[1], values[2] );
                         angleDom[idom] = OrientedAngle(p0, p1, gref, g); // -pi<angle<+pi
                         //MESSAGE("  angle=" << angleDom[idom]);
                       }
@@ -11552,7 +11745,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   std::map<int, std::map<long,int> > nodeQuadDomains;
   std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
 
-  MESSAGE(".. Creation of elements: simple junction");
+  //MESSAGE(".. Creation of elements: simple junction");
   if (createJointElems)
   {
     int idg;
@@ -11603,7 +11796,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   //     iterate on mutipleNodesToFace
   //     iterate on edgesMultiDomains
 
-  MESSAGE(".. Creation of elements: multiple junction");
+  //MESSAGE(".. Creation of elements: multiple junction");
   if (createJointElems)
   {
     // --- iterate on mutipleNodesToFace
@@ -11676,7 +11869,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   faceOrEdgeDom.clear();
   feDom.clear();
 
-  MESSAGE(".. Modification of elements");
+  //MESSAGE(".. Modification of elements");
   for (int idomain = idom0; idomain < nbDomains; idomain++)
   {
     std::map<int, std::map<int, int> >::const_iterator itnod = nodeDomains.begin();
@@ -11702,7 +11895,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
               feDom[vtkId] = idomain;
               faceOrEdgeDom[aCell] = emptyMap;
               faceOrEdgeDom[aCell][idomain] = vtkId; // affect face or edge to the first domain only
-              //MESSAGE("affect cell " << this->GetMeshDS()->fromVtkToSmds(vtkId) << " domain " << idomain
+              //MESSAGE("affect cell " << this->GetMeshDS()->FromVtkToSmds(vtkId) << " domain " << idomain
               //        << " type " << vtkType << " downId " << downId);
             }
       }
@@ -11734,7 +11927,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       {
         int idom = itdom->first;
         int vtkVolId = itdom->second;
-        //MESSAGE("modify nodes of cell " << this->GetMeshDS()->fromVtkToSmds(vtkVolId) << " domain " << idom);
+        //MESSAGE("modify nodes of cell " << this->GetMeshDS()->FromVtkToSmds(vtkVolId) << " domain " << idom);
         localClonedNodeIds.clear();
         for (itn = oldNodes.begin(); itn != oldNodes.end(); ++itn)
         {
@@ -11759,7 +11952,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   }
 
   meshDS->CleanDownWardConnectivity(); // Mesh has been modified, downward connectivity is no more usable, free memory
-  grid->BuildLinks();
+  grid->DeleteLinks();
 
   CHRONOSTOP(50);
   counters::stats();
@@ -11778,9 +11971,9 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
  */
 bool SMESH_MeshEditor::CreateFlatElementsOnFacesGroups(const std::vector<TIDSortedElemSet>& theElems)
 {
-  MESSAGE("-------------------------------------------------");
-  MESSAGE("SMESH_MeshEditor::CreateFlatElementsOnFacesGroups");
-  MESSAGE("-------------------------------------------------");
+  // MESSAGE("-------------------------------------------------");
+  // MESSAGE("SMESH_MeshEditor::CreateFlatElementsOnFacesGroups");
+  // MESSAGE("-------------------------------------------------");
 
   SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS();
 
@@ -11801,8 +11994,7 @@ bool SMESH_MeshEditor::CreateFlatElementsOnFacesGroups(const std::vector<TIDSort
     TIDSortedElemSet::const_iterator elemItr = domain.begin();
     for ( ; elemItr != domain.end(); ++elemItr )
     {
-      SMDS_MeshElement* anElem = (SMDS_MeshElement*) *elemItr;
-      SMDS_MeshFace* aFace = dynamic_cast<SMDS_MeshFace*> (anElem);
+      const SMDS_MeshFace* aFace = meshDS->DownCast<SMDS_MeshFace> ( *elemItr );
       if (!aFace)
         continue;
       // MESSAGE("aFace=" << aFace->GetID());
@@ -11811,11 +12003,11 @@ bool SMESH_MeshEditor::CreateFlatElementsOnFacesGroups(const std::vector<TIDSort
 
       // --- clone the nodes, create intermediate nodes for non medium nodes of a quad face
 
-      SMDS_ElemIteratorPtr nodeIt = aFace->nodesIterator();
+      SMDS_NodeIteratorPtr nodeIt = aFace->nodeIterator();
       while (nodeIt->more())
       {
-        const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*> (nodeIt->next());
-        bool isMedium = isQuad && (aFace->IsMediumNode(node));
+        const SMDS_MeshNode* node = nodeIt->next();
+        bool isMedium = ( isQuad && aFace->IsMediumNode( node ));
         if (isMedium)
           ln2.push_back(node);
         else
@@ -11923,7 +12115,7 @@ bool SMESH_MeshEditor::CreateFlatElementsOnFacesGroups(const std::vector<TIDSort
 
       // --- modify the face
 
-      aFace->ChangeNodes(&ln[0], ln.size());
+      const_cast<SMDS_MeshFace*>( aFace )->ChangeNodes( &ln[0], ln.size() );
     }
   }
   return true;
@@ -11942,9 +12134,9 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
                                       std::vector<double>&            nodesCoords,
                                       std::vector<std::vector<int> >& listOfListOfNodes)
 {
-  MESSAGE("--------------------------------");
-  MESSAGE("SMESH_MeshEditor::CreateHoleSkin");
-  MESSAGE("--------------------------------");
+  // MESSAGE("--------------------------------");
+  // MESSAGE("SMESH_MeshEditor::CreateHoleSkin");
+  // MESSAGE("--------------------------------");
 
   // --- zone of volumes to remove is given :
   //     1 either by a geom shape (one or more vertices) and a radius,
@@ -12034,7 +12226,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
 
   SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS();
   meshDS->BuildDownWardConnectivity(true);
-  SMDS_UnstructuredGrid* grid = meshDS->getGrid();
+  SMDS_UnstructuredGrid* grid = meshDS->GetGrid();
 
   // --- set of volumes detected inside
 
@@ -12046,7 +12238,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
 
   if (isNodeGroup) // --- a group of nodes is provided : find all the volumes using one or more of this nodes
   {
-    MESSAGE("group of nodes provided");
+    //MESSAGE("group of nodes provided");
     SMDS_ElemIteratorPtr elemIt = groupDS->GetElements();
     while ( elemIt->more() )
     {
@@ -12061,14 +12253,14 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
       while (volItr->more())
       {
         vol = (SMDS_MeshElement*)volItr->next();
-        setOfInsideVol.insert(vol->getVtkId());
+        setOfInsideVol.insert(vol->GetVtkID());
         sgrp->Add(vol->GetID());
       }
     }
   }
   else if (isNodeCoords)
   {
-    MESSAGE("list of nodes coordinates provided");
+    //MESSAGE("list of nodes coordinates provided");
     size_t i = 0;
     int k = 0;
     while ( i < nodesCoords.size()-2 )
@@ -12078,13 +12270,13 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
       double z = nodesCoords[i++];
       gp_Pnt p = gp_Pnt(x, y ,z);
       gpnts.push_back(p);
-      MESSAGE("TopoDS_Vertex " << k << " " << p.X() << " " << p.Y() << " " << p.Z());
+      //MESSAGE("TopoDS_Vertex " << k << " " << p.X() << " " << p.Y() << " " << p.Z());
       k++;
     }
   }
   else // --- no group, no coordinates : use the vertices of the geom shape provided, and radius
   {
-    MESSAGE("no group of nodes provided, using vertices from geom shape, and radius");
+    //MESSAGE("no group of nodes provided, using vertices from geom shape, and radius");
     TopTools_IndexedMapOfShape vertexMap;
     TopExp::MapShapes( theShape, TopAbs_VERTEX, vertexMap );
     gp_Pnt p = gp_Pnt(0,0,0);
@@ -12096,20 +12288,17 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
       const TopoDS_Vertex& vertex = TopoDS::Vertex( vertexMap( i ));
       p = BRep_Tool::Pnt(vertex);
       gpnts.push_back(p);
-      MESSAGE("TopoDS_Vertex " << i << " " << p.X() << " " << p.Y() << " " << p.Z());
+      //MESSAGE("TopoDS_Vertex " << i << " " << p.X() << " " << p.Y() << " " << p.Z());
     }
   }
 
   if (gpnts.size() > 0)
   {
-    int nodeId = 0;
     const SMDS_MeshNode* startNode = theNodeSearcher->FindClosestTo(gpnts[0]);
-    if (startNode)
-      nodeId = startNode->GetID();
-    MESSAGE("nodeId " << nodeId);
+    //MESSAGE("startNode->nodeId " << nodeId);
 
     double radius2 = radius*radius;
-    MESSAGE("radius2 " << radius2);
+    //MESSAGE("radius2 " << radius2);
 
     // --- volumes on start node
 
@@ -12119,7 +12308,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
     while (volItr->more())
     {
       startVol = (SMDS_MeshElement*)volItr->next();
-      setOfVolToCheck.insert(startVol->getVtkId());
+      setOfVolToCheck.insert(startVol->GetVtkID());
     }
     if (setOfVolToCheck.empty())
     {
@@ -12138,7 +12327,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
     {
       std::set<int>::iterator it = setOfVolToCheck.begin();
       int vtkId = *it;
-      MESSAGE("volume to check,  vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+      //MESSAGE("volume to check,  vtkId " << vtkId << " smdsId " << meshDS->FromVtkToSmds(vtkId));
       bool volInside = false;
       vtkIdType npts = 0;
       vtkIdType* pts = 0;
@@ -12149,7 +12338,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
         if (mapOfNodeDistance2.count(pts[i]))
         {
           distance2 = mapOfNodeDistance2[pts[i]];
-          MESSAGE("point " << pts[i] << " distance2 " << distance2);
+          //MESSAGE("point " << pts[i] << " distance2 " << distance2);
         }
         else
         {
@@ -12167,19 +12356,19 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
             }
           }
           mapOfNodeDistance2[pts[i]] = distance2;
-          MESSAGE("  point "  << pts[i]  << " distance2 " << distance2 << " coords " << coords[0] << " " << coords[1] << " " <<  coords[2]);
+          //MESSAGE("  point "  << pts[i]  << " distance2 " << distance2 << " coords " << coords[0] << " " << coords[1] << " " <<  coords[2]);
         }
         if (distance2 < radius2)
         {
           volInside = true; // one or more nodes inside the domain
-          sgrp->Add(meshDS->fromVtkToSmds(vtkId));
+          sgrp->Add(meshDS->FromVtkToSmds(vtkId));
           break;
         }
       }
       if (volInside)
       {
         setOfInsideVol.insert(vtkId);
-        MESSAGE("  volume inside,  vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+        //MESSAGE("  volume inside,  vtkId " << vtkId << " smdsId " << meshDS->FromVtkToSmds(vtkId));
         int neighborsVtkIds[NBMAXNEIGHBORS];
         int downIds[NBMAXNEIGHBORS];
         unsigned char downTypes[NBMAXNEIGHBORS];
@@ -12191,7 +12380,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
       else
       {
         setOfOutsideVol.insert(vtkId);
-        MESSAGE("  volume outside, vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+        //MESSAGE("  volume outside, vtkId " << vtkId << " smdsId " << meshDS->FromVtkToSmds(vtkId));
       }
       setOfVolToCheck.erase(vtkId);
     }
@@ -12204,7 +12393,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
   std::set<int> setOfVolToReCheck;
   while (addedInside)
   {
-    MESSAGE(" --------------------------- re check");
+    //MESSAGE(" --------------------------- re check");
     addedInside = false;
     std::set<int>::iterator itv = setOfInsideVol.begin();
     for (; itv != setOfInsideVol.end(); ++itv)
@@ -12226,7 +12415,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
       int vtkId = *it;
       if (grid->GetCellType(vtkId) == VTK_HEXAHEDRON)
       {
-        MESSAGE("volume to recheck,  vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+        //MESSAGE("volume to recheck,  vtkId " << vtkId << " smdsId " << meshDS->FromVtkToSmds(vtkId));
         int countInside = 0;
         int neighborsVtkIds[NBMAXNEIGHBORS];
         int downIds[NBMAXNEIGHBORS];
@@ -12235,12 +12424,12 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
         for (int n = 0; n < nbNeighbors; n++)
           if (setOfInsideVol.count(neighborsVtkIds[n]))
             countInside++;
-        MESSAGE("countInside " << countInside);
+        //MESSAGE("countInside " << countInside);
         if (countInside > 1)
         {
-          MESSAGE("  volume inside,  vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+          //MESSAGE("  volume inside,  vtkId " << vtkId << " smdsId " << meshDS->FromVtkToSmds(vtkId));
           setOfInsideVol.insert(vtkId);
-          sgrp->Add(meshDS->fromVtkToSmds(vtkId));
+          sgrp->Add(meshDS->FromVtkToSmds(vtkId));
           addedInside = true;
         }
         else
@@ -12261,7 +12450,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
   for (; it != setOfInsideVol.end(); ++it)
   {
     int vtkId = *it;
-    //MESSAGE("  vtkId " << vtkId  << " smdsId " << meshDS->fromVtkToSmds(vtkId));
+    //MESSAGE("  vtkId " << vtkId  << " smdsId " << meshDS->FromVtkToSmds(vtkId));
     int neighborsVtkIds[NBMAXNEIGHBORS];
     int downIds[NBMAXNEIGHBORS];
     unsigned char downTypes[NBMAXNEIGHBORS];
@@ -12280,7 +12469,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
         int vtkFaceId = grid->getDownArray(downTypes[n])->getVtkCellId(downIds[n]);
         if (vtkFaceId >= 0)
         {
-          sgrpi->Add(meshDS->fromVtkToSmds(vtkFaceId));
+          sgrpi->Add(meshDS->FromVtkToSmds(vtkFaceId));
           // find also the smds edges on this face
           int nbEdges = grid->getDownArray(downTypes[n])->getNumberOfDownCells(downIds[n]);
           const int* dEdges = grid->getDownArray(downTypes[n])->getDownCells(downIds[n]);
@@ -12289,7 +12478,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
           {
             int vtkEdgeId = grid->getDownArray(dTypes[i])->getVtkCellId(dEdges[i]);
             if (vtkEdgeId >= 0)
-              sgrpei->Add(meshDS->fromVtkToSmds(vtkEdgeId));
+              sgrpei->Add(meshDS->FromVtkToSmds(vtkEdgeId));
           }
         }
       }
@@ -12299,7 +12488,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
         skinFaces[face] = vtkId;
         int vtkFaceId = grid->getDownArray(downTypes[n])->getVtkCellId(downIds[n]);
         if (vtkFaceId >= 0)
-          sgrps->Add(meshDS->fromVtkToSmds(vtkFaceId));
+          sgrps->Add(meshDS->FromVtkToSmds(vtkFaceId));
       }
     }
   }
@@ -12318,7 +12507,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
   {
     const SMDS_MeshElement *elem = itelem->next();
     int shapeId = elem->getshapeId();
-    int vtkId = elem->getVtkId();
+    int   vtkId = elem->GetVtkID();
     if (!shapeIdToVtkIdSet.count(shapeId))
     {
       shapeIdToVtkIdSet[shapeId] = emptySet;
@@ -12335,7 +12524,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
   for (; itShape != shapeIdToVtkIdSet.end(); ++itShape)
   {
     int shapeId = itShape->first;
-    MESSAGE(" --- Shape ID --- "<< shapeId);
+    //MESSAGE(" --- Shape ID --- "<< shapeId);
     shapeIdToEdges[shapeId] = emptyEdges;
 
     std::vector<int> nodesEdges;
@@ -12344,7 +12533,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
     for (; its != itShape->second.end(); ++its)
     {
       int vtkId = *its;
-      MESSAGE("     " << vtkId);
+      //MESSAGE("     " << vtkId);
       int neighborsVtkIds[NBMAXNEIGHBORS];
       int downIds[NBMAXNEIGHBORS];
       unsigned char downTypes[NBMAXNEIGHBORS];
@@ -12353,7 +12542,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
       {
         if (neighborsVtkIds[n]<0) // only smds faces are considered as neighbors here
           continue;
-        int smdsId = meshDS->fromVtkToSmds(neighborsVtkIds[n]);
+        int smdsId = meshDS->FromVtkToSmds(neighborsVtkIds[n]);
         const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
         if ( shapeIds.count(elem->getshapeId()) && !sgrps->Contains(elem)) // edge : neighbor in the set of shape, not in the group
         {
@@ -12365,7 +12554,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
             int nbNodes = grid->getDownArray(downTypes[n])->getNodes(downIds[n],vtkNodeId);
             nodesEdges.push_back(vtkNodeId[0]);
             nodesEdges.push_back(vtkNodeId[nbNodes-1]);
-            MESSAGE("       --- nodes " << vtkNodeId[0]+1 << " " << vtkNodeId[nbNodes-1]+1);
+            //MESSAGE("       --- nodes " << vtkNodeId[0]+1 << " " << vtkNodeId[nbNodes-1]+1);
           }
         }
       }
@@ -12375,9 +12564,9 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
     order.clear();
     if (nodesEdges.size() > 0)
     {
-      order.push_back(nodesEdges[0]); MESSAGE("       --- back " << order.back()+1); // SMDS id = VTK id + 1;
+      order.push_back(nodesEdges[0]); //MESSAGE("       --- back " << order.back()+1); // SMDS id = VTK id + 1;
       nodesEdges[0] = -1;
-      order.push_back(nodesEdges[1]); MESSAGE("       --- back " << order.back()+1);
+      order.push_back(nodesEdges[1]); //MESSAGE("       --- back " << order.back()+1);
       nodesEdges[1] = -1; // do not reuse this edge
       bool found = true;
       while (found)
@@ -12396,7 +12585,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
               found = false;
             else
             {
-              order.push_back(nodesEdges[i-1]); MESSAGE("       --- back " << order.back()+1);
+              order.push_back(nodesEdges[i-1]); //MESSAGE("       --- back " << order.back()+1);
               nodesEdges[i-1] = -1;
             }
           else // even ==> use the next one
@@ -12404,7 +12593,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
               found = false;
             else
             {
-              order.push_back(nodesEdges[i+1]); MESSAGE("       --- back " << order.back()+1);
+              order.push_back(nodesEdges[i+1]); //MESSAGE("       --- back " << order.back()+1);
               nodesEdges[i+1] = -1;
             }
         }
@@ -12426,7 +12615,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
             found = false;
           else
           {
-            order.push_front(nodesEdges[i-1]); MESSAGE("       --- front " << order.front()+1);
+            order.push_front(nodesEdges[i-1]); //MESSAGE("       --- front " << order.front()+1);
             nodesEdges[i-1] = -1;
           }
         else // even ==> use the next one
@@ -12434,7 +12623,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
             found = false;
           else
           {
-            order.push_front(nodesEdges[i+1]); MESSAGE("       --- front " << order.front()+1);
+            order.push_front(nodesEdges[i+1]); //MESSAGE("       --- front " << order.front()+1);
             nodesEdges[i+1] = -1;
           }
       }
@@ -12447,7 +12636,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
     for (; itl != order.end(); itl++)
     {
       nodes.push_back((*itl) + 1); // SMDS id = VTK id + 1;
-      MESSAGE("              ordered node " << nodes[nodes.size()-1]);
+      //MESSAGE("              ordered node " << nodes[nodes.size()-1]);
     }
     listOfListOfNodes.push_back(nodes);
   }
@@ -12507,7 +12696,7 @@ bool SMESH_MeshEditor::Make2DMeshFrom3D()
       // add new face based on volume nodes
       if (aMesh->FindElement( nodes, SMDSAbs_Face, /*noMedium=*/false) )
       {
-        nbExisted++; // face already exsist
+        nbExisted++; // face already exists
       }
       else
       {
@@ -12591,9 +12780,9 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
 
   SMDS_ElemIteratorPtr eIt;
   if (elements.empty()) eIt = aMesh->elementsIterator(elemType);
-  else                  eIt = elemSetIterator( elements );
+  else                  eIt = SMESHUtils::elemSetIterator( elements );
 
-  while (eIt->more())
+  while ( eIt->more() )
   {
     const SMDS_MeshElement* elem = eIt->next();
     const int              iQuad = elem->IsQuadratic();
@@ -12674,7 +12863,7 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
     else if ( elem->GetType() == SMDSAbs_Face ) // elem is a face ------------------------
     {
       avoidSet.clear(), avoidSet.insert( elem );
-      elemNodes.assign( SMDS_MeshElement::iterator( elem->interlacedNodesElemIterator() ),
+      elemNodes.assign( SMDS_MeshElement::iterator( elem->interlacedNodesIterator() ),
                         SMDS_MeshElement::iterator() );
       elemNodes.push_back( elemNodes[0] );
       nodes.resize( 2 + iQuad );
@@ -12707,7 +12896,7 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
         tgtNodes.resize( srcNodes.size() );
         for ( inode = 0; inode < srcNodes.size(); ++inode )
           tgtNodes[inode] = getNodeWithSameID( tgtMeshDS, srcNodes[inode] );
-        if ( aroundElements && tgtEditor.GetMeshDS()->FindElement( tgtNodes,
+        if ( /*aroundElements && */tgtEditor.GetMeshDS()->FindElement( tgtNodes,
                                                                    missType,
                                                                    /*noMedium=*/false))
           continue;
@@ -12718,7 +12907,7 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
       for ( size_t i = 0; i < missingBndElems.size(); ++i )
       {
         TConnectivity& nodes = missingBndElems[ i ];
-        if ( aroundElements && tgtEditor.GetMeshDS()->FindElement( nodes,
+        if ( /*aroundElements && */tgtEditor.GetMeshDS()->FindElement( nodes,
                                                                    missType,
                                                                    /*noMedium=*/false))
           continue;
@@ -12770,7 +12959,7 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
     else // store present elements to add them to a group
       for ( size_t i = 0 ; i < presentBndElems.size(); ++i )
       {
-        presentEditor->myLastCreatedElems.Append( presentBndElems[ i ]);
+        presentEditor->myLastCreatedElems.push_back( presentBndElems[ i ]);
       }
 
   } // loop on given elements
@@ -12781,11 +12970,11 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
   if ( group )
   {
     if ( SMESHDS_Group* g = dynamic_cast<SMESHDS_Group*>( group->GetGroupDS() ))
-      for ( int i = 0; i < tgtEditor.myLastCreatedElems.Size(); ++i )
-        g->SMDSGroup().Add( tgtEditor.myLastCreatedElems( i+1 ));
+      for ( size_t i = 0; i < tgtEditor.myLastCreatedElems.size(); ++i )
+        g->SMDSGroup().Add( tgtEditor.myLastCreatedElems[ i ]);
   }
-  tgtEditor.myLastCreatedElems.Clear();
-  tgtEditor2.myLastCreatedElems.Clear();
+  tgtEditor.myLastCreatedElems.clear();
+  tgtEditor2.myLastCreatedElems.clear();
 
   // -----------------------
   // 5. Copy given elements
@@ -12793,7 +12982,7 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
   if ( toCopyElements && targetMesh != myMesh )
   {
     if (elements.empty()) eIt = aMesh->elementsIterator(elemType);
-    else                  eIt = elemSetIterator( elements );
+    else                  eIt = SMESHUtils::elemSetIterator( elements );
     while (eIt->more())
     {
       const SMDS_MeshElement* elem = eIt->next();
@@ -12802,7 +12991,7 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
         tgtNodes[inode] = getNodeWithSameID( tgtMeshDS, elem->GetNode(inode) );
       tgtEditor.AddElement( tgtNodes, elemToCopy.Init( elem ));
 
-      tgtEditor.myLastCreatedElems.Clear();
+      tgtEditor.myLastCreatedElems.clear();
     }
   }
   return nbAddedBnd;
@@ -12828,7 +13017,7 @@ void SMESH_MeshEditor::copyPosition( const SMDS_MeshNode* from,
 
   case SMDS_TOP_FACE:
   {
-    const SMDS_FacePosition* fPos = static_cast< const SMDS_FacePosition* >( pos );
+    SMDS_FacePositionPtr fPos = pos;
     GetMeshDS()->SetNodeOnFace( to, from->getshapeId(),
                                 fPos->GetUParameter(), fPos->GetVParameter() );
     break;
@@ -12836,7 +13025,7 @@ void SMESH_MeshEditor::copyPosition( const SMDS_MeshNode* from,
   case SMDS_TOP_EDGE:
   {
     // WARNING: it is dangerous to set equal nodes on one EDGE!!!!!!!!
-    const SMDS_EdgePosition* ePos = static_cast< const SMDS_EdgePosition* >( pos );
+    SMDS_EdgePositionPtr ePos = pos;
     GetMeshDS()->SetNodeOnEdge( to, from->getshapeId(), ePos->GetUParameter() );
     break;
   }
@@ -12849,3 +13038,482 @@ void SMESH_MeshEditor::copyPosition( const SMDS_MeshNode* from,
   default:;
   }
 }
+
+namespace // utils for MakePolyLine
+{
+  //================================================================================
+  /*!
+   * \brief Sequence of found points and a current point data
+   */
+  struct Path
+  {
+    std::vector< gp_XYZ >   myPoints;
+    double                  myLength;
+
+    int                     mySrcPntInd; //!< start point index
+    const SMDS_MeshElement* myFace;
+    SMESH_NodeXYZ           myNode1;
+    SMESH_NodeXYZ           myNode2;
+    int                     myNodeInd1;
+    int                     myNodeInd2;
+    double                  myDot1;
+    double                  myDot2;
+    TIDSortedElemSet        myElemSet, myAvoidSet;
+
+    Path(): myLength(0.0), myFace(0) {}
+
+    bool SetCutAtCorner( const SMESH_NodeXYZ&    cornerNode,
+                         const SMDS_MeshElement* face,
+                         const gp_XYZ&           plnNorm,
+                         const gp_XYZ&           plnOrig );
+
+    void AddPoint( const gp_XYZ& p );
+
+    bool Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig );
+
+    bool ReachSamePoint( const Path& other );
+
+    static void Remove( std::vector< Path > & paths, size_t& i );
+  };
+
+  //================================================================================
+  /*!
+   * \brief Return true if this Path meats another
+   */
+  //================================================================================
+
+  bool Path::ReachSamePoint( const Path& other )
+  {
+    return ( mySrcPntInd != other.mySrcPntInd &&
+             myFace == other.myFace );
+  }
+
+  //================================================================================
+  /*!
+   * \brief Remove a path from a vector
+   */
+  //================================================================================
+
+  void Path::Remove( std::vector< Path > & paths, size_t& i )
+  {
+    if ( paths.size() > 1 )
+    {
+      size_t j = paths.size() - 1; // last item to be removed
+      if ( i < j )
+      {
+        paths[ i ].myPoints.swap( paths[ j ].myPoints );
+        paths[ i ].myLength    = paths[ j ].myLength;
+        paths[ i ].mySrcPntInd = paths[ j ].mySrcPntInd;
+        paths[ i ].myFace      = paths[ j ].myFace;
+        paths[ i ].myNode1     = paths[ j ].myNode1;
+        paths[ i ].myNode2     = paths[ j ].myNode2;
+        paths[ i ].myNodeInd1  = paths[ j ].myNodeInd1;
+        paths[ i ].myNodeInd2  = paths[ j ].myNodeInd2;
+        paths[ i ].myDot1      = paths[ j ].myDot1;
+        paths[ i ].myDot2      = paths[ j ].myDot2;
+      }
+    }
+    paths.pop_back();
+    if ( i > 0 )
+      --i;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Store a point that is at a node of a face if the face is intersected by plane.
+   *        Return false if the node is a sole intersection point of the face and the plane
+   */
+  //================================================================================
+
+  bool Path::SetCutAtCorner( const SMESH_NodeXYZ&    cornerNode,
+                             const SMDS_MeshElement* face,
+                             const gp_XYZ&           plnNorm,
+                             const gp_XYZ&           plnOrig )
+  {
+    if ( face == myFace )
+      return false;
+    myNodeInd1 = face->GetNodeIndex( cornerNode._node );
+    myNodeInd2 = ( myNodeInd1 + 1 ) % face->NbCornerNodes();
+    int   ind3 = ( myNodeInd1 + 2 ) % face->NbCornerNodes();
+    myNode1.Set( face->GetNode( ind3 ));
+    myNode2.Set( face->GetNode( myNodeInd2 ));
+
+    myDot1 = plnNorm * ( myNode1 - plnOrig );
+    myDot2 = plnNorm * ( myNode2 - plnOrig );
+
+    bool ok = ( myDot1 * myDot2 < 0 );
+    if ( !ok && myDot1 * myDot2 == 0 )
+    {
+      ok = ( myDot1 != myDot2 );
+      if ( ok && myFace )
+        ok = ( myFace->GetNodeIndex(( myDot1 == 0 ? myNode1 : myNode2 )._node ) < 0 );
+    }
+    if ( ok )
+    {
+      myFace = face;
+      myDot1 = 0;
+      AddPoint( cornerNode );
+    }
+    return ok;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Store a point and update myLength
+   */
+  //================================================================================
+
+  void Path::AddPoint( const gp_XYZ& p )
+  {
+    if ( !myPoints.empty() )
+      myLength += ( p - myPoints.back() ).Modulus();
+    else
+      myLength = 0;
+    myPoints.push_back( p );
+  }
+
+  //================================================================================
+  /*!
+   * \brief Try to find the next point
+   *  \param [in] plnNorm - cutting plane normal
+   *  \param [in] plnOrig - cutting plane origin
+   */
+  //================================================================================
+
+  bool Path::Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig )
+  {
+    int nodeInd3 = ( myNodeInd1 + 1 ) % myFace->NbCornerNodes();
+    if ( myNodeInd2 == nodeInd3 )
+      nodeInd3 = ( myNodeInd1 + 2 ) % myFace->NbCornerNodes();
+
+    SMESH_NodeXYZ node3 = myFace->GetNode( nodeInd3 );
+    double         dot3 = plnNorm * ( node3 - plnOrig );
+
+    if ( dot3 * myDot1 < 0. )
+    {
+      myNode2    = node3;
+      myNodeInd2 = nodeInd3;
+      myDot2     = dot3;
+    }
+    else if ( dot3 * myDot2 < 0. )
+    {
+      myNode1    = node3;
+      myNodeInd1 = nodeInd3;
+      myDot1     = dot3;
+    }
+    else if ( dot3 == 0. )
+    {
+      SMDS_ElemIteratorPtr fIt = node3._node->GetInverseElementIterator(SMDSAbs_Face);
+      while ( fIt->more() )
+        if ( SetCutAtCorner( node3, fIt->next(), plnNorm, plnOrig ))
+          return true;
+      return false;
+    }
+    else if ( myDot2 == 0. )
+    {
+      SMESH_NodeXYZ node2 = myNode2; // copy as myNode2 changes in SetCutAtCorner()
+      SMDS_ElemIteratorPtr fIt = node2._node->GetInverseElementIterator(SMDSAbs_Face);
+      while ( fIt->more() )
+        if ( SetCutAtCorner( node2, fIt->next(), plnNorm, plnOrig ))
+          return true;
+      return false;
+    }
+
+    double r = Abs( myDot1 / ( myDot2 - myDot1 ));
+    AddPoint( myNode1 * ( 1 - r ) + myNode2 * r );
+
+    myAvoidSet.clear();
+    myAvoidSet.insert( myFace );
+    myFace = SMESH_MeshAlgos::FindFaceInSet( myNode1._node, myNode2._node,
+                                             myElemSet,   myAvoidSet,
+                                             &myNodeInd1, &myNodeInd2 );
+    return myFace;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Compute a path between two points of PolySegment
+   */
+  struct PolyPathCompute
+  {
+    SMESH_MeshEditor::TListOfPolySegments& mySegments; //!< inout PolySegment's
+    std::vector< Path >&                   myPaths;    //!< path of each of segments to compute
+    SMESH_Mesh*                            myMesh;
+    mutable std::vector< std::string >     myErrors;
+
+    PolyPathCompute( SMESH_MeshEditor::TListOfPolySegments& theSegments,
+                     std::vector< Path >&                   thePaths,
+                     SMESH_Mesh*                            theMesh):
+      mySegments( theSegments ),
+      myPaths( thePaths ),
+      myMesh( theMesh ),
+      myErrors( theSegments.size() )
+    {
+    }
+#undef SMESH_CAUGHT
+#define SMESH_CAUGHT myErrors[i] =
+    void operator() ( const int i ) const
+    {
+      SMESH_TRY;
+      const_cast< PolyPathCompute* >( this )->Compute( i );
+      SMESH_CATCH( SMESH::returnError );
+    }
+#undef SMESH_CAUGHT
+    //================================================================================
+    /*!
+     * \brief Compute a path of a given segment
+     */
+    //================================================================================
+
+    void Compute( const int iSeg )
+    {
+      SMESH_MeshEditor::PolySegment& polySeg = mySegments[ iSeg ];
+
+      // get a cutting plane
+
+      gp_XYZ p1 = SMESH_NodeXYZ( polySeg.myNode1[0] );
+      gp_XYZ p2 = SMESH_NodeXYZ( polySeg.myNode1[1] );
+      if ( polySeg.myNode2[0] ) p1 = 0.5 * ( p1 + SMESH_NodeXYZ( polySeg.myNode2[0] ));
+      if ( polySeg.myNode2[1] ) p2 = 0.5 * ( p2 + SMESH_NodeXYZ( polySeg.myNode2[1] ));
+
+      gp_XYZ plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
+      gp_XYZ plnOrig = p2;
+
+      // find paths connecting the 2 end points of polySeg
+
+      std::vector< Path > paths; paths.reserve(10);
+
+      // initialize paths
+
+      for ( int iP = 0; iP < 2; ++iP ) // loop on the polySeg end points
+      {
+        Path path;
+        path.mySrcPntInd = iP;
+        size_t nbPaths = paths.size();
+
+        if ( polySeg.myNode2[ iP ] && polySeg.myNode2[ iP ] != polySeg.myNode1[ iP ] )
+        {
+          while (( path.myFace = SMESH_MeshAlgos::FindFaceInSet( polySeg.myNode1[ iP ],
+                                                                 polySeg.myNode2[ iP ],
+                                                                 path.myElemSet,
+                                                                 path.myAvoidSet,
+                                                                 &path.myNodeInd1,
+                                                                 &path.myNodeInd2 )))
+          {
+            path.myNode1.Set( polySeg.myNode1[ iP ]);
+            path.myNode2.Set( polySeg.myNode2[ iP ]);
+            path.myDot1 = plnNorm * ( path.myNode1 - plnOrig );
+            path.myDot2 = plnNorm * ( path.myNode2 - plnOrig );
+            path.myPoints.clear();
+            path.AddPoint( 0.5 * ( path.myNode1 + path.myNode2 ));
+            path.myAvoidSet.insert( path.myFace );
+            paths.push_back( path );
+          }
+          if ( nbPaths == paths.size() )
+            throw SALOME_Exception ( SMESH_Comment("No face edge found by point ") << iP+1
+                                     << " in a PolySegment " << iSeg );
+        }
+        else // an end point is at node
+        {
+          std::set<const SMDS_MeshNode* > nodes;
+          SMDS_ElemIteratorPtr fIt = polySeg.myNode1[ iP ]->GetInverseElementIterator(SMDSAbs_Face);
+          while ( fIt->more() )
+          {
+            path.myPoints.clear();
+            if ( path.SetCutAtCorner( polySeg.myNode1[ iP ], fIt->next(), plnNorm, plnOrig ))
+            {
+              if (( path.myDot1 * path.myDot2 != 0 ) ||
+                  ( nodes.insert( path.myDot1 == 0 ? path.myNode1._node : path.myNode2._node ).second ))
+                paths.push_back( path );
+            }
+          }
+        }
+
+        // look for a one-segment path
+        for ( size_t i = 0; i < nbPaths; ++i )
+          for ( size_t j = nbPaths; j < paths.size(); ++j )
+            if ( paths[i].myFace == paths[j].myFace )
+            {
+              myPaths[ iSeg ].myPoints.push_back( paths[i].myPoints[0] );
+              myPaths[ iSeg ].myPoints.push_back( paths[j].myPoints[0] );
+              paths.clear();
+            }
+      }
+
+      // extend paths
+
+      myPaths[ iSeg ].myLength = 1e100;
+
+      while ( paths.size() >= 2 )
+      {
+        for ( size_t i = 0; i < paths.size(); ++i )
+        {
+          Path& path = paths[ i ];
+          if ( !path.Extend( plnNorm, plnOrig ) ||         // path reached a mesh boundary
+               path.myLength > myPaths[ iSeg ].myLength ) // path is longer than others
+          {
+            Path::Remove( paths, i );
+            continue;
+          }
+
+          // join paths that reach same point
+          for ( size_t j = 0; j < paths.size(); ++j )
+          {
+            if ( i != j && paths[i].ReachSamePoint( paths[j] ))
+            {
+              double   distLast = ( paths[i].myPoints.back() - paths[j].myPoints.back() ).Modulus();
+              double fullLength = ( paths[i].myLength + paths[j].myLength + distLast );
+              if ( fullLength < myPaths[ iSeg ].myLength )
+              {
+                myPaths[ iSeg ].myLength = fullLength;
+                std::vector< gp_XYZ > & allPoints = myPaths[ iSeg ].myPoints;
+                allPoints.swap( paths[i].myPoints );
+                allPoints.insert( allPoints.end(),
+                                  paths[j].myPoints.rbegin(),
+                                  paths[j].myPoints.rend() );
+              }
+              Path::Remove( paths, i );
+              Path::Remove( paths, j );
+            }
+          }
+        }
+        if ( !paths.empty() && (int) paths[0].myPoints.size() > myMesh->NbFaces() )
+          throw SALOME_Exception(LOCALIZED( "Infinite loop in MakePolyLine()"));
+      }
+
+      if ( myPaths[ iSeg ].myPoints.empty() )
+        throw SALOME_Exception( SMESH_Comment("Can't find a full path for PolySegment #") << iSeg );
+
+    } // PolyPathCompute::Compute()
+
+  }; // struct PolyPathCompute
+
+} // namespace
+
+//=======================================================================
+//function : MakePolyLine
+//purpose  : Create a polyline consisting of 1D mesh elements each lying on a 2D element of
+//           the initial mesh
+//=======================================================================
+
+void SMESH_MeshEditor::MakePolyLine( TListOfPolySegments&   theSegments,
+                                     SMESHDS_Group*         theGroup,
+                                     SMESH_ElementSearcher* theSearcher)
+{
+  std::vector< Path > segPaths( theSegments.size() ); // path of each of segments
+
+  SMESH_ElementSearcher* searcher = theSearcher;
+  SMESHUtils::Deleter<SMESH_ElementSearcher> delSearcher;
+  if ( !searcher )
+  {
+    searcher = SMESH_MeshAlgos::GetElementSearcher( *GetMeshDS() );
+    delSearcher._obj = searcher;
+  }
+
+  // get cutting planes
+
+  std::vector< bool > isVectorOK( theSegments.size(), true );
+  const double planarCoef = 0.333; // plane height in planar case
+
+  for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
+  {
+    PolySegment& polySeg = theSegments[ iSeg ];
+
+    gp_XYZ p1 = SMESH_NodeXYZ( polySeg.myNode1[0] );
+    gp_XYZ p2 = SMESH_NodeXYZ( polySeg.myNode1[1] );
+    if ( polySeg.myNode2[0] ) p1 = 0.5 * ( p1 + SMESH_NodeXYZ( polySeg.myNode2[0] ));
+    if ( polySeg.myNode2[1] ) p2 = 0.5 * ( p2 + SMESH_NodeXYZ( polySeg.myNode2[1] ));
+
+    gp_XYZ plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
+
+    isVectorOK[ iSeg ] = ( plnNorm.Modulus() > std::numeric_limits<double>::min() );
+    if ( !isVectorOK[ iSeg ])
+    {
+      gp_XYZ pMid = 0.5 * ( p1 + p2 );
+      const SMDS_MeshElement* face;
+      polySeg.myMidProjPoint = searcher->Project( pMid, SMDSAbs_Face, &face );
+      polySeg.myVector       = polySeg.myMidProjPoint.XYZ() - pMid;
+
+      gp_XYZ faceNorm;
+      SMESH_MeshAlgos::FaceNormal( face, faceNorm );
+
+      if ( polySeg.myVector.Magnitude() < Precision::Confusion() ||
+           polySeg.myVector * faceNorm  < Precision::Confusion() )
+      {
+        polySeg.myVector = faceNorm;
+        polySeg.myMidProjPoint = pMid + faceNorm * ( p1 - p2 ).Modulus() * planarCoef;
+      }
+    }
+    else
+    {
+      polySeg.myVector = plnNorm ^ ( p1 - p2 );
+    }
+  }
+
+  // assure that inverse elements are constructed, avoid their concurrent building in threads
+  GetMeshDS()->nodesIterator()->next()->NbInverseElements();
+
+  // find paths
+
+  PolyPathCompute algo( theSegments, segPaths, myMesh );
+  OSD_Parallel::For( 0, theSegments.size(), algo, theSegments.size() == 1 );
+
+  for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
+    if ( !algo.myErrors[ iSeg ].empty() )
+      throw SALOME_Exception( algo.myErrors[ iSeg ].c_str() );
+
+  // create an 1D mesh
+
+  const SMDS_MeshNode *n, *nPrev = 0;
+  SMESHDS_Mesh* mesh = GetMeshDS();
+
+  for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
+  {
+    const Path& path = segPaths[iSeg];
+    if ( path.myPoints.size() < 2 )
+      continue;
+
+    double tol = path.myLength / path.myPoints.size() / 1000.;
+    if ( !nPrev || ( SMESH_NodeXYZ( nPrev ) - path.myPoints[0] ).SquareModulus() > tol*tol )
+    {
+      nPrev = mesh->AddNode( path.myPoints[0].X(), path.myPoints[0].Y(), path.myPoints[0].Z() );
+      myLastCreatedNodes.push_back( nPrev );
+    }
+    for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
+    {
+      n = mesh->AddNode( path.myPoints[iP].X(), path.myPoints[iP].Y(), path.myPoints[iP].Z() );
+      myLastCreatedNodes.push_back( n );
+
+      const SMDS_MeshElement* elem = mesh->AddEdge( nPrev, n );
+      myLastCreatedElems.push_back( elem );
+      if ( theGroup )
+        theGroup->Add( elem );
+
+      nPrev = n;
+    }
+
+    // return a vector
+
+    gp_XYZ pMid = 0.5 * ( path.myPoints[0] + path.myPoints.back() );
+    if ( isVectorOK[ iSeg ])
+    {
+      // find the most distance point of a path
+      double maxDist = 0;
+      for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
+      {
+        double dist = Abs( theSegments[iSeg].myVector * ( path.myPoints[iP] - path.myPoints[0] ));
+        if ( dist > maxDist )
+        {
+          maxDist = dist;
+          theSegments[iSeg].myMidProjPoint = path.myPoints[iP];
+        }
+      }
+      if ( maxDist < Precision::Confusion() ) // planar case
+        theSegments[iSeg].myMidProjPoint =
+          pMid + theSegments[iSeg].myVector.XYZ().Normalized() * path.myLength * planarCoef;
+    }
+    theSegments[iSeg].myVector = gp_Vec( pMid, theSegments[iSeg].myMidProjPoint );
+  }
+
+  return;
+}