Salome HOME
untabify
[modules/smesh.git] / src / SMESH / SMESH_MeshEditor.cxx
index b858c372f1249ad88d97ff39fb0c5659a2cc07d2..e84487512375142c0e588951c39b798ffab2bd3b 100644 (file)
@@ -1,6 +1,6 @@
-//  SMESH SMESH : idl implementation based on 'SMESH' unit's classes
+//  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
 //
-//  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+//  Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
 //
 //  This library is free software; you can redistribute it and/or
 //  License along with this library; if not, write to the Free Software
 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
-// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-//
-//
+//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
+//  SMESH SMESH : idl implementation based on 'SMESH' unit's classes
 // File      : SMESH_MeshEditor.cxx
 // Created   : Mon Apr 12 16:10:22 2004
 // Author    : Edward AGAPOV (eap)
-
-
+//
 #include "SMESH_MeshEditor.hxx"
 
 #include "SMDS_FaceOfNodes.hxx"
@@ -49,6 +47,7 @@
 #include "utilities.h"
 
 #include <BRep_Tool.hxx>
+#include <BRepClass3d_SolidClassifier.hxx>
 #include <ElCLib.hxx>
 #include <Extrema_GenExtPS.hxx>
 #include <Extrema_POnSurf.hxx>
 #include <GeomAdaptor_Surface.hxx>
 #include <Geom_Curve.hxx>
 #include <Geom_Surface.hxx>
+#include <Precision.hxx>
 #include <TColStd_ListOfInteger.hxx>
+#include <TopAbs_State.hxx>
 #include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
 #include <TopTools_ListIteratorOfListOfShape.hxx>
 #include <TopTools_ListOfShape.hxx>
+#include <TopTools_SequenceOfShape.hxx>
 #include <TopoDS.hxx>
 #include <TopoDS_Face.hxx>
 #include <gp.hxx>
@@ -76,6 +78,8 @@
 
 #include <map>
 #include <set>
+#include <numeric>
+#include <limits>
 
 #define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
 
@@ -88,24 +92,23 @@ typedef map<const SMDS_MeshElement*, list<const SMDS_MeshElement*> > TElemOfElem
 //typedef TNodeOfNodeVecMap::iterator                                  TNodeOfNodeVecMapItr;
 //typedef map<const SMDS_MeshElement*, vector<TNodeOfNodeVecMapItr> >  TElemOfVecOfMapNodesMap;
 
-struct TNodeXYZ : public gp_XYZ {
-  TNodeXYZ( const SMDS_MeshNode* n ):gp_XYZ( n->X(), n->Y(), n->Z() ) {}
-};
-
-typedef pair< const SMDS_MeshNode*, const SMDS_MeshNode* > NLink;
-
 //=======================================================================
 /*!
- * \brief A sorted pair of nodes
+ * \brief SMDS_MeshNode -> gp_XYZ convertor
  */
 //=======================================================================
 
-struct TLink: public NLink
+struct TNodeXYZ : public gp_XYZ
 {
-  TLink(const SMDS_MeshNode* n1, const SMDS_MeshNode* n2 ):NLink( n1, n2 )
-  { if ( n1->GetID() < n2->GetID() ) std::swap( first, second ); }
-  TLink(const NLink& link ):NLink( link )
-  { if ( first->GetID() < second->GetID() ) std::swap( first, second ); }
+  TNodeXYZ( const SMDS_MeshNode* n ):gp_XYZ( n->X(), n->Y(), n->Z() ) {}
+  double Distance( const SMDS_MeshNode* n )
+  {
+    return gp_Vec( *this, TNodeXYZ( n )).Magnitude();
+  }
+  double SquareDistance( const SMDS_MeshNode* n )
+  {
+    return gp_Vec( *this, TNodeXYZ( n )).SquareMagnitude();
+  }
 };
 
 //=======================================================================
@@ -283,17 +286,17 @@ bool SMESH_MeshEditor::Remove (const list< int >& theIDs,
             smmap.insert( sm );
     }
     // Find sub-meshes to notify about modification
-//     SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
-//     while ( nodeIt->more() ) {
-//       const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
-//       const SMDS_PositionPtr& aPosition = node->GetPosition();
-//       if ( aPosition.get() ) {
-//         if ( int aShapeID = aPosition->GetShapeId() ) {
-//           if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
-//             smmap.insert( sm );
-//         }
-//       }
-//     }
+    //     SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
+    //     while ( nodeIt->more() ) {
+    //       const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+    //       const SMDS_PositionPtr& aPosition = node->GetPosition();
+    //       if ( aPosition.get() ) {
+    //         if ( int aShapeID = aPosition->GetShapeId() ) {
+    //           if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
+    //             smmap.insert( sm );
+    //         }
+    //       }
+    //     }
 
     // Do remove
     if ( isNodes )
@@ -309,9 +312,9 @@ bool SMESH_MeshEditor::Remove (const list< int >& theIDs,
       (*smIt)->ComputeStateEngine( SMESH_subMesh::MESH_ENTITY_REMOVED );
   }
 
-//   // Check if the whole mesh becomes empty
-//   if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( 1 ) )
-//     sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
+  //   // Check if the whole mesh becomes empty
+  //   if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( 1 ) )
+  //     sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
 
   return true;
 }
@@ -1171,11 +1174,13 @@ void SMESH_MeshEditor::AddToSameGroups (const SMDS_MeshElement* elemToAdd,
                                         SMESHDS_Mesh *          aMesh)
 {
   const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
-  set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
-  for ( ; grIt != groups.end(); grIt++ ) {
-    SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
-    if ( group && group->Contains( elemInGroups ))
-      group->SMDSGroup().Add( elemToAdd );
+  if (!groups.empty()) {
+    set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
+    for ( ; grIt != groups.end(); grIt++ ) {
+      SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
+      if ( group && group->Contains( elemInGroups ))
+        group->SMDSGroup().Add( elemToAdd );
+    }
   }
 }
 
@@ -1200,6 +1205,25 @@ void SMESH_MeshEditor::RemoveElemFromGroups (const SMDS_MeshElement* removeelem,
   }
 }
 
+//=======================================================================
+//function : ReplaceElemInGroups
+//purpose  : replace elemToRm by elemToAdd in the all groups
+//=======================================================================
+
+void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement* elemToRm,
+                                            const SMDS_MeshElement* elemToAdd,
+                                            SMESHDS_Mesh *          aMesh)
+{
+  const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
+  if (!groups.empty()) {
+    set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
+    for ( ; grIt != groups.end(); grIt++ ) {
+      SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
+      if ( group && group->SMDSGroup().Remove( elemToRm ) && elemToAdd )
+        group->SMDSGroup().Add( elemToAdd );
+    }
+  }
+}
 
 //=======================================================================
 //function : QuadToTri
@@ -1404,7 +1428,7 @@ double getAngle(const SMDS_MeshElement * tr1,
 // and able to return nodes by that ID
 // =================================================
 class LinkID_Gen {
- public:
+public:
 
   LinkID_Gen( const SMESHDS_Mesh* theMesh )
     :myMesh( theMesh ), myMaxID( theMesh->MaxNodeID() + 1)
@@ -1427,7 +1451,7 @@ class LinkID_Gen {
     return true;
   }
 
- private:
+private:
   LinkID_Gen();
   const SMESHDS_Mesh* myMesh;
   long                myMaxID;
@@ -1460,10 +1484,10 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
   // 1. map of elements with their linkIDs
   // 2. map of linkIDs with their elements
 
-  map< TLink, list< const SMDS_MeshElement* > > mapLi_listEl;
-  map< TLink, list< const SMDS_MeshElement* > >::iterator itLE;
-  map< const SMDS_MeshElement*, set< TLink > >  mapEl_setLi;
-  map< const SMDS_MeshElement*, set< TLink > >::iterator itEL;
+  map< SMESH_TLink, list< const SMDS_MeshElement* > > mapLi_listEl;
+  map< SMESH_TLink, list< const SMDS_MeshElement* > >::iterator itLE;
+  map< const SMDS_MeshElement*, set< SMESH_TLink > >  mapEl_setLi;
+  map< const SMDS_MeshElement*, set< SMESH_TLink > >::iterator itEL;
 
   TIDSortedElemSet::iterator itElem;
   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
@@ -1482,7 +1506,7 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
 
     // fill maps
     for ( i = 0; i < 3; i++ ) {
-      TLink link( aNodes[i], aNodes[i+1] );
+      SMESH_TLink link( aNodes[i], aNodes[i+1] );
       // check if elements sharing a link can be fused
       itLE = mapLi_listEl.find( link );
       if ( itLE != mapLi_listEl.end() ) {
@@ -1508,7 +1532,7 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
     int nbElems = (*itLE).second.size();
     if ( nbElems < 2  ) {
       const SMDS_MeshElement* elem = (*itLE).second.front();
-      TLink link = (*itLE).first;
+      SMESH_TLink link = (*itLE).first;
       mapEl_setLi[ elem ].erase( link );
       if ( mapEl_setLi[ elem ].empty() )
         mapEl_setLi.erase( elem );
@@ -1534,11 +1558,11 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
 
     // search elements to fuse starting from startElem or links of elements
     // fused earlyer - startLinks
-    list< TLink > startLinks;
+    list< SMESH_TLink > startLinks;
     while ( startElem || !startLinks.empty() ) {
       while ( !startElem && !startLinks.empty() ) {
         // Get an element to start, by a link
-        TLink linkId = startLinks.front();
+        SMESH_TLink linkId = startLinks.front();
         startLinks.pop_front();
         itLE = mapLi_listEl.find( linkId );
         if ( itLE != mapLi_listEl.end() ) {
@@ -1554,15 +1578,15 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
       if ( startElem ) {
         // Get candidates to be fused
         const SMDS_MeshElement *tr1 = startElem, *tr2 = 0, *tr3 = 0;
-        const TLink *link12, *link13;
+        const SMESH_TLink *link12, *link13;
         startElem = 0;
         ASSERT( mapEl_setLi.find( tr1 ) != mapEl_setLi.end() );
-        set< TLink >& setLi = mapEl_setLi[ tr1 ];
+        set< SMESH_TLink >& setLi = mapEl_setLi[ tr1 ];
         ASSERT( !setLi.empty() );
-        set< TLink >::iterator itLi;
+        set< SMESH_TLink >::iterator itLi;
         for ( itLi = setLi.begin(); itLi != setLi.end(); itLi++ )
         {
-          const TLink & link = (*itLi);
+          const SMESH_TLink & link = (*itLi);
           itLE = mapLi_listEl.find( link );
           if ( itLE == mapLi_listEl.end() )
             continue;
@@ -1583,10 +1607,10 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
           }
 
           // add other links of elem to list of links to re-start from
-          set< TLink >& links = mapEl_setLi[ elem ];
-          set< TLink >::iterator it;
+          set< SMESH_TLink >& links = mapEl_setLi[ elem ];
+          set< SMESH_TLink >::iterator it;
           for ( it = links.begin(); it != links.end(); it++ ) {
-            const TLink& link2 = (*it);
+            const SMESH_TLink& link2 = (*it);
             if ( link2 != link )
               startLinks.push_back( link2 );
           }
@@ -1736,15 +1760,15 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
 //=============================================================================
 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 ] << ")");
+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 ] << ")");
 }
 
 //=======================================================================
@@ -1755,7 +1779,7 @@ static void swap( int i1, int i2, int idNodes[], gp_Pnt P[] )
 //=======================================================================
 
 int SMESH_MeshEditor::SortQuadNodes (const SMDS_Mesh * theMesh,
-                                     int               idNodes[] )
+int               idNodes[] )
 {
   gp_Pnt P[4];
   int i;
@@ -1784,10 +1808,10 @@ int SMESH_MeshEditor::SortQuadNodes (const SMDS_Mesh * theMesh,
       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());
-//     }
+    //     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;
 }
@@ -1903,7 +1927,7 @@ bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh,
           faceNodes.insert( idNodes[ 2 ] );
           faceNodes.insert( idNodes[ iMin ] );
           DUMPSO( "loop " << iLoop2 << " id2 " << idNodes[ 1 ] << " id3 " << idNodes[ 2 ]
-            << " leastDist = " << leastDist);
+                  << " leastDist = " << leastDist);
           if ( leastDist <= DBL_MIN )
             break;
         }
@@ -1941,11 +1965,11 @@ bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh,
     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());
-//    }
+  //   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.;
@@ -1997,11 +2021,11 @@ bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh,
     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]);
-//   }
+  //   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;
 }*/
@@ -2009,11 +2033,11 @@ bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh,
 //================================================================================
 /*!
  * \brief Return nodes linked to the given one
 * \param theNode - the node
 * \param linkedNodes - the found nodes
 * \param type - the type of elements to check
 *
 * Medium nodes are ignored
+ * \param theNode - the node
+ * \param linkedNodes - the found nodes
+ * \param type - the type of elements to check
+ *
+ * Medium nodes are ignored
  */
 //================================================================================
 
@@ -2047,8 +2071,8 @@ void SMESH_MeshEditor::GetLinkedNodes( const SMDS_MeshNode* theNode,
             iAfter  = SMESH_MesherHelper::WrapIndex( iAfter, nb );
             iBefore = SMESH_MesherHelper::WrapIndex( iBefore, nb );
           }
-          linkedNodes.insert( elem->GetNode( iAfter ));
-          linkedNodes.insert( elem->GetNode( iBefore ));
+          linkedNodes.insert( elem->GetNodeWrap( iAfter ));
+          linkedNodes.insert( elem->GetNodeWrap( iBefore ));
         }
       }
     }
@@ -2293,14 +2317,14 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
     }
     int nbElemOnFace = 0;
     itElem = theElems.begin();
-     // loop on not yet smoothed elements: look for elems on a face
+    // loop on not yet smoothed elements: look for elems on a face
     while ( itElem != theElems.end() ) {
       if ( faceSubMesh && nbElemOnFace == faceSubMesh->NbElements() )
         break; // all elements found
 
       const SMDS_MeshElement* elem = *itElem;
       if ( !elem || elem->GetType() != SMDSAbs_Face || elem->NbNodes() < 3 ||
-          ( faceSubMesh && !faceSubMesh->Contains( elem ))) {
+           ( faceSubMesh && !faceSubMesh->Contains( elem ))) {
         ++itElem;
         continue;
       }
@@ -2359,19 +2383,19 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
         if ( uvMap.find( node ) == uvMap.end() )
           uvCheckNodes.push_back( node );
         // add nodes of elems sharing node
-//         SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Face);
-//         while ( eIt->more() ) {
-//           const SMDS_MeshElement* e = eIt->next();
-//           if ( e != elem ) {
-//             SMDS_ElemIteratorPtr nIt = e->nodesIterator();
-//             while ( nIt->more() ) {
-//               const SMDS_MeshNode* n =
-//                 static_cast<const SMDS_MeshNode*>( nIt->next() );
-//               if ( uvMap.find( n ) == uvMap.end() )
-//                 uvCheckNodes.push_back( n );
-//             }
-//           }
-//         }
+        //         SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Face);
+        //         while ( eIt->more() ) {
+        //           const SMDS_MeshElement* e = eIt->next();
+        //           if ( e != elem ) {
+        //             SMDS_ElemIteratorPtr nIt = e->nodesIterator();
+        //             while ( nIt->more() ) {
+        //               const SMDS_MeshNode* n =
+        //                 static_cast<const SMDS_MeshNode*>( nIt->next() );
+        //               if ( uvMap.find( n ) == uvMap.end() )
+        //                 uvCheckNodes.push_back( n );
+        //             }
+        //           }
+        //         }
       }
       // check UV on face
       list< const SMDS_MeshNode* >::iterator n = uvCheckNodes.begin();
@@ -2443,9 +2467,8 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
     // fix nodes on mesh boundary
 
     if ( checkBoundaryNodes ) {
-      typedef pair<const SMDS_MeshNode*, const SMDS_MeshNode*> TLink;
-      map< TLink, int > linkNbMap; // how many times a link encounters in elemsOnFace
-      map< TLink, int >::iterator link_nb;
+      map< NLink, int > linkNbMap; // how many times a link encounters in elemsOnFace
+      map< NLink, int >::iterator link_nb;
       // put all elements links to linkNbMap
       list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
       for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
@@ -2454,10 +2477,10 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
         if(elem->IsQuadratic())
           nbn = nbn/2;
         // loop on elem links: insert them in linkNbMap
-        const SMDS_MeshNode* curNode, *prevNode = elem->GetNode( nbn );
+        const SMDS_MeshNode* curNode, *prevNode = elem->GetNodeWrap( nbn );
         for ( int iN = 0; iN < nbn; ++iN ) {
           curNode = elem->GetNode( iN );
-          TLink link;
+          NLink link;
           if ( curNode < prevNode ) link = make_pair( curNode , prevNode );
           else                      link = make_pair( prevNode , curNode );
           prevNode = curNode;
@@ -2772,60 +2795,57 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
 
   int iNode, nbSame = 0, iNotSameNode = 0, iSameNode = 0;
   vector<int> sames(nbNodes);
-
-  //bool issimple[nbNodes];
   vector<bool> issimple(nbNodes);
 
   for ( iNode = 0; iNode < nbNodes; iNode++ ) {
     TNodeOfNodeListMapItr nnIt = newNodesItVec[ iNode ];
     const SMDS_MeshNode*                 node         = nnIt->first;
     const list< const SMDS_MeshNode* > & listNewNodes = nnIt->second;
-    if ( listNewNodes.empty() )
+    if ( listNewNodes.empty() ) {
       return;
-
-    if(listNewNodes.size()==nbSteps) {
-      issimple[iNode] = true;
-    }
-    else {
-      issimple[iNode] = false;
     }
 
+    issimple[iNode] = (listNewNodes.size()==nbSteps);
+
     itNN[ iNode ] = listNewNodes.begin();
     prevNod[ iNode ] = node;
     nextNod[ iNode ] = listNewNodes.front();
-//cout<<"iNode="<<iNode<<endl;
-//cout<<" prevNod[iNode]="<< prevNod[iNode]<<" nextNod[iNode]="<< nextNod[iNode]<<endl;
-    if ( prevNod[ iNode ] != nextNod [ iNode ])
-      iNotSameNode = iNode;
-    else {
-      iSameNode = iNode;
-      //nbSame++;
-      sames[nbSame++] = iNode;
+    if( !issimple[iNode] ) {
+      if ( prevNod[ iNode ] != nextNod [ iNode ])
+        iNotSameNode = iNode;
+      else {
+        iSameNode = iNode;
+        //nbSame++;
+        sames[nbSame++] = iNode;
+      }
     }
   }
-//cout<<"1 nbSame="<<nbSame<<endl;
+
+  //cout<<"  nbSame = "<<nbSame<<endl;
   if ( nbSame == nbNodes || nbSame > 2) {
-    MESSAGE( " Too many same nodes of element " << elem->GetID() );
+    //MESSAGE( " Too many same nodes of element " << elem->GetID() );
+    INFOS( " Too many same nodes of element " << elem->GetID() );
     return;
   }
 
-//  if( elem->IsQuadratic() && nbSame>0 ) {
-//    MESSAGE( "Can not rotate quadratic element " << elem->GetID() );
-//    return;
-//  }
+  //  if( elem->IsQuadratic() && nbSame>0 ) {
+  //    MESSAGE( "Can not rotate quadratic element " << elem->GetID() );
+  //    return;
+  //  }
 
   int iBeforeSame = 0, iAfterSame = 0, iOpposSame = 0;
+  int nbBaseNodes = ( elem->IsQuadratic() ? nbNodes/2 : nbNodes );
   if ( nbSame > 0 ) {
-    iBeforeSame = ( iSameNode == 0 ? nbNodes - 1 : iSameNode - 1 );
-    iAfterSame  = ( iSameNode + 1 == nbNodes ? 0 : iSameNode + 1 );
+    iBeforeSame = ( iSameNode == 0 ? nbBaseNodes - 1 : iSameNode - 1 );
+    iAfterSame  = ( iSameNode + 1 == nbBaseNodes ? 0 : iSameNode + 1 );
     iOpposSame  = ( iSameNode - 2 < 0  ? iSameNode + 2 : iSameNode - 2 );
   }
 
-//if(nbNodes==8)
-//cout<<" prevNod[0]="<< prevNod[0]<<" prevNod[1]="<< prevNod[1]
-//    <<" prevNod[2]="<< prevNod[2]<<" prevNod[3]="<< prevNod[4]
-//    <<" prevNod[4]="<< prevNod[4]<<" prevNod[5]="<< prevNod[5]
-//    <<" prevNod[6]="<< prevNod[6]<<" prevNod[7]="<< prevNod[7]<<endl;
+  //if(nbNodes==8)
+  //cout<<" prevNod[0]="<< prevNod[0]<<" prevNod[1]="<< prevNod[1]
+  //    <<" prevNod[2]="<< prevNod[2]<<" prevNod[3]="<< prevNod[4]
+  //    <<" prevNod[4]="<< prevNod[4]<<" prevNod[5]="<< prevNod[5]
+  //    <<" prevNod[6]="<< prevNod[6]<<" prevNod[7]="<< prevNod[7]<<endl;
 
   // check element orientation
   int i0 = 0, i2 = 2;
@@ -2833,11 +2853,8 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
     //MESSAGE("Reversed elem " << elem );
     i0 = 2;
     i2 = 0;
-    if ( nbSame > 0 ) {
-      int iAB = iAfterSame + iBeforeSame;
-      iBeforeSame = iAB - iBeforeSame;
-      iAfterSame  = iAB - iAfterSame;
-    }
+    if ( nbSame > 0 )
+      std::swap( iBeforeSame, iAfterSame );
   }
 
   // make new elements
@@ -2856,10 +2873,9 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
           nextNod[ iNode ] = *itNN[ iNode ];
           itNN[ iNode ]++;
         }
-        else if(!elem->IsQuadratic() ||
-           elem->IsQuadratic() && elem->IsMediumNode(prevNod[iNode]) ) {
+        else if(!elem->IsQuadratic() || elem->IsMediumNode(prevNod[iNode]) ) {
           // we have to use each second node
-          itNN[ iNode ]++;
+          //itNN[ iNode ]++;
           nextNod[ iNode ] = *itNN[ iNode ];
           itNN[ iNode ]++;
         }
@@ -2918,8 +2934,9 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
                                       midlNod[0], nextNod[2], midlNod[1], prevNod[2]);
           }
           else if(nbSame==1) { // quadratic triangle
-            if(sames[0]==2)
+            if(sames[0]==2) {
               return; // medium node on axis
+            }
             else if(sames[0]==0) {
               aNewElem = aMesh->AddFace(prevNod[0], nextNod[1], prevNod[1],
                                         nextNod[2], midlNod[1], prevNod[2]);
@@ -2929,8 +2946,9 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
                                         midlNod[0], nextNod[2], prevNod[2]);
             }
           }
-          else
+          else {
             return;
+          }
         }
         break;
       }
@@ -2965,37 +2983,139 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
       }
       case 6: { // quadratic triangle
         // create pentahedron with 15 nodes
-        if(i0>0) { // reversed case
-          aNewElem = aMesh->AddVolume (prevNod[0], prevNod[2], prevNod[1],
-                                       nextNod[0], nextNod[2], nextNod[1],
-                                       prevNod[5], prevNod[4], prevNod[3],
-                                       nextNod[5], nextNod[4], nextNod[3],
-                                       midlNod[0], midlNod[2], midlNod[1]);
-        }
-        else { // not reversed case
-          aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2],
-                                       nextNod[0], nextNod[1], nextNod[2],
-                                       prevNod[3], prevNod[4], prevNod[5],
-                                       nextNod[3], nextNod[4], nextNod[5],
-                                       midlNod[0], midlNod[1], midlNod[2]);
+        if(nbSame==0) {
+          if(i0>0) { // reversed case
+            aNewElem = aMesh->AddVolume (prevNod[0], prevNod[2], prevNod[1],
+                                         nextNod[0], nextNod[2], nextNod[1],
+                                         prevNod[5], prevNod[4], prevNod[3],
+                                         nextNod[5], nextNod[4], nextNod[3],
+                                         midlNod[0], midlNod[2], midlNod[1]);
+          }
+          else { // not reversed case
+            aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2],
+                                         nextNod[0], nextNod[1], nextNod[2],
+                                         prevNod[3], prevNod[4], prevNod[5],
+                                         nextNod[3], nextNod[4], nextNod[5],
+                                         midlNod[0], midlNod[1], midlNod[2]);
+          }
+        }
+        else if(nbSame==1) {
+          // 2d order pyramid of 13 nodes
+          //SMDS_MeshVolume* AddVolumeWithID(int n1, int n2, int n3, int n4, int n5,
+          //                                 int n12,int n23,int n34,int n41,
+          //                                 int n15,int n25,int n35,int n45, int ID);
+          int n5 = iSameNode;
+          int n1,n4,n41,n15,n45;
+          if(i0>0) { // reversed case
+            n1 = ( n5 + 1 == nbBaseNodes ? 0 : n5 + 1 );
+            n4 = ( n5 == 0 ? nbBaseNodes - 1 : n5 - 1 );
+            n41 = n1 + 3;
+            n15 = n5 + 3;
+            n45 = n4 + 3;
+          }
+          else {
+            n1 = ( n5 == 0 ? nbBaseNodes - 1 : n5 - 1 );
+            n4 = ( n5 + 1 == nbBaseNodes ? 0 : n5 + 1 );
+            n41 = n4 + 3;
+            n15 = n1 + 3;
+            n45 = n5 + 3;
+          }
+          aNewElem = aMesh->AddVolume(prevNod[n1], nextNod[n1],
+                                      nextNod[n4], prevNod[n4], prevNod[n5],
+                                      midlNod[n1], nextNod[n41],
+                                      midlNod[n4], prevNod[n41],
+                                      prevNod[n15], nextNod[n15],
+                                      nextNod[n45], prevNod[n45]);
+        }
+        else if(nbSame==2) {
+          // 2d order tetrahedron of 10 nodes
+          //SMDS_MeshVolume* AddVolumeWithID(int n1, int n2, int n3, int n4,
+          //                                 int n12,int n23,int n31,
+          //                                 int n14,int n24,int n34, int ID);
+          int n1 = iNotSameNode;
+          int n2,n3,n12,n23,n31;
+          if(i0>0) { // reversed case
+            n2 = ( n1 == 0 ? nbBaseNodes - 1 : n1 - 1 );
+            n3 = ( n1 + 1 == nbBaseNodes ? 0 : n1 + 1 );
+            n12 = n2 + 3;
+            n23 = n3 + 3;
+            n31 = n1 + 3;
+          }
+          else {
+            n2 = ( n1 + 1 == nbBaseNodes ? 0 : n1 + 1 );
+            n3 = ( n1 == 0 ? nbBaseNodes - 1 : n1 - 1 );
+            n12 = n1 + 3;
+            n23 = n2 + 3;
+            n31 = n3 + 3;
+          }
+          aNewElem = aMesh->AddVolume (prevNod[n1], prevNod[n2], prevNod[n3], nextNod[n1],
+                                       prevNod[n12], prevNod[n23], prevNod[n31],
+                                       midlNod[n1], nextNod[n12], nextNod[n31]);
         }
         break;
       }
       case 8: { // quadratic quadrangle
-        // create hexahedron with 20 nodes
-        if(i0>0) { // reversed case
-          aNewElem = aMesh->AddVolume (prevNod[0], prevNod[3], prevNod[2], prevNod[1],
-                                       nextNod[0], nextNod[3], nextNod[2], nextNod[1],
-                                       prevNod[7], prevNod[6], prevNod[5], prevNod[4],
-                                       nextNod[7], nextNod[6], nextNod[5], nextNod[4],
-                                       midlNod[0], midlNod[3], midlNod[2], midlNod[1]);
-        }
-        else { // not reversed case
-          aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2], prevNod[3],
-                                       nextNod[0], nextNod[1], nextNod[2], nextNod[3],
-                                       prevNod[4], prevNod[5], prevNod[6], prevNod[7],
-                                       nextNod[4], nextNod[5], nextNod[6], nextNod[7],
-                                       midlNod[0], midlNod[1], midlNod[2], midlNod[3]);
+        if(nbSame==0) {
+          // create hexahedron with 20 nodes
+          if(i0>0) { // reversed case
+            aNewElem = aMesh->AddVolume (prevNod[0], prevNod[3], prevNod[2], prevNod[1],
+                                         nextNod[0], nextNod[3], nextNod[2], nextNod[1],
+                                         prevNod[7], prevNod[6], prevNod[5], prevNod[4],
+                                         nextNod[7], nextNod[6], nextNod[5], nextNod[4],
+                                         midlNod[0], midlNod[3], midlNod[2], midlNod[1]);
+          }
+          else { // not reversed case
+            aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2], prevNod[3],
+                                         nextNod[0], nextNod[1], nextNod[2], nextNod[3],
+                                         prevNod[4], prevNod[5], prevNod[6], prevNod[7],
+                                         nextNod[4], nextNod[5], nextNod[6], nextNod[7],
+                                         midlNod[0], midlNod[1], midlNod[2], midlNod[3]);
+          }
+        }
+        else if(nbSame==1) { 
+          // --- pyramid + pentahedron - can not be created since it is needed 
+          // additional middle node ot the center of face
+          INFOS( " Sweep for face " << elem->GetID() << " can not be created" );
+          return;
+        }
+        else if(nbSame==2) {
+          // 2d order Pentahedron with 15 nodes
+          //SMDS_MeshVolume* AddVolumeWithID(int n1, int n2, int n3, int n4, int n5, int n6,
+          //                                 int n12,int n23,int n31,int n45,int n56,int n64,
+          //                                 int n14,int n25,int n36, int ID);
+          int n1,n2,n4,n5;
+          if ( prevNod[ iBeforeSame ] == nextNod[ iBeforeSame ] ) {
+            // iBeforeSame is same too
+            n1 = iBeforeSame;
+            n2 = iOpposSame;
+            n4 = iSameNode;
+            n5 = iAfterSame;
+          }
+          else {
+            // iAfterSame is same too
+            n1 = iSameNode;
+            n2 = iBeforeSame;
+            n4 = iAfterSame;
+            n5 = iOpposSame;
+          }
+          int n12,n45,n14,n25;
+          if(i0>0) { //reversed case
+            n12 = n1 + 4;
+            n45 = n5 + 4;
+            n14 = n4 + 4;
+            n25 = n2 + 4;
+          }
+          else {
+            n12 = n2 + 4;
+            n45 = n4 + 4;
+            n14 = n1 + 4;
+            n25 = n5 + 4;
+          }
+          aNewElem = aMesh->AddVolume (prevNod[n1], prevNod[n2], nextNod[n2],
+                                       prevNod[n4], prevNod[n5], nextNod[n5],
+                                       prevNod[n12], midlNod[n2], nextNod[n12],
+                                       prevNod[n45], midlNod[n5], nextNod[n45],
+                                       prevNod[n14], prevNod[n25], nextNod[n25]);
         }
         break;
       }
@@ -3285,7 +3405,7 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
               const SMDS_MeshFace * f = aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ]);
               if ( !f )
                 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ));
-              else if ( nodes[ 1 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
+              else if ( nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
                 aMesh->ChangeElementNodes( f, nodes, nbn );
               break;
             }
@@ -3293,7 +3413,7 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
               const SMDS_MeshFace * f = aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]);
               if ( !f )
                 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ));
-              else if ( nodes[ 1 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
+              else if ( nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
                 aMesh->ChangeElementNodes( f, nodes, nbn );
               break;
             }
@@ -3302,20 +3422,40 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
                 if(nbn==6) { /////// quadratic triangle
                   const SMDS_MeshFace * f = aMesh->FindFace( nodes[0], nodes[2], nodes[4],
                                                              nodes[1], nodes[3], nodes[5] );
-                  if ( !f )
+                  if ( !f ) {
                     myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4],
                                                              nodes[1], nodes[3], nodes[5]));
-                  else if ( nodes[ 2 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
-                    aMesh->ChangeElementNodes( f, nodes, nbn );
+                  }
+                  else if ( nodes[ 2 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 )) {
+                    const SMDS_MeshNode** tmpnodes = new const SMDS_MeshNode*[6];
+                    tmpnodes[0] = nodes[0];
+                    tmpnodes[1] = nodes[2];
+                    tmpnodes[2] = nodes[4];
+                    tmpnodes[3] = nodes[1];
+                    tmpnodes[4] = nodes[3];
+                    tmpnodes[5] = nodes[5];
+                    aMesh->ChangeElementNodes( f, tmpnodes, nbn );
+                  }
                 }
                 else {       /////// quadratic quadrangle
                   const SMDS_MeshFace * f = aMesh->FindFace( nodes[0], nodes[2], nodes[4], nodes[6],
                                                              nodes[1], nodes[3], nodes[5], nodes[7] );
-                  if ( !f )
+                  if ( !f ) {
                     myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
                                                              nodes[1], nodes[3], nodes[5], nodes[7]));
-                  else if ( nodes[ 2 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
-                    aMesh->ChangeElementNodes( f, nodes, nbn );
+                  }
+                  else if ( nodes[ 2 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 )) {
+                    const SMDS_MeshNode** tmpnodes = new const SMDS_MeshNode*[8];
+                    tmpnodes[0] = nodes[0];
+                    tmpnodes[1] = nodes[2];
+                    tmpnodes[2] = nodes[4];
+                    tmpnodes[3] = nodes[6];
+                    tmpnodes[4] = nodes[1];
+                    tmpnodes[5] = nodes[3];
+                    tmpnodes[6] = nodes[5];
+                    tmpnodes[7] = nodes[7];
+                    aMesh->ChangeElementNodes( f, tmpnodes, nbn );
+                  }
                 }
               }
               else { //////// polygon
@@ -3323,7 +3463,7 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
                 const SMDS_MeshFace * f = aMesh->FindFace( polygon_nodes );
                 if ( !f )
                   myLastCreatedElems.Append(aMesh->AddPolygonalFace(polygon_nodes));
-                else if ( nodes[ 1 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
+                else if ( nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
                   aMesh->ChangeElementNodes( f, nodes, nbn );
               }
             }
@@ -3436,20 +3576,25 @@ SMESH_MeshEditor::RotationSweep(TIDSortedElemSet & theElems,
 
     // loop on elem nodes
     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
-    while ( itN->more() )
-    {
+    while ( itN->more() ) {
       // check if a node has been already sweeped
       const SMDS_MeshNode* node = cast2Node( itN->next() );
+
+      gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
+      double coord[3];
+      aXYZ.Coord( coord[0], coord[1], coord[2] );
+      bool isOnAxis = ( aLine.SquareDistance( aXYZ ) <= aSqTol );
+
       TNodeOfNodeListMapItr nIt = mapNewNodes.find( node );
       if ( nIt == mapNewNodes.end() ) {
         nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
         list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
 
         // make new nodes
-        gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
-        double coord[3];
-        aXYZ.Coord( coord[0], coord[1], coord[2] );
-        bool isOnAxis = ( aLine.SquareDistance( aXYZ ) <= aSqTol );
+        //gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
+        //double coord[3];
+        //aXYZ.Coord( coord[0], coord[1], coord[2] );
+        //bool isOnAxis = ( aLine.SquareDistance( aXYZ ) <= aSqTol );
         const SMDS_MeshNode * newNode = node;
         for ( int i = 0; i < theNbSteps; i++ ) {
           if ( !isOnAxis ) {
@@ -3470,37 +3615,52 @@ SMESH_MeshEditor::RotationSweep(TIDSortedElemSet & theElems,
             newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
             myLastCreatedNodes.Append(newNode);
             srcNodes.Append( node );
+            listNewNodes.push_back( newNode );
+          }
+          else {
+            listNewNodes.push_back( newNode );
+            if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
+              listNewNodes.push_back( newNode );
+            }
           }
-          listNewNodes.push_back( newNode );
         }
       }
-      else {
+      /*
+        else {
         // if current elem is quadratic and current node is not medium
         // we have to check - may be it is needed to insert additional nodes
         if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
-          list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
-          if(listNewNodes.size()==theNbSteps) {
-            listNewNodes.clear();
-            // make new nodes
-            gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
-            double coord[3];
-            aXYZ.Coord( coord[0], coord[1], coord[2] );
-            const SMDS_MeshNode * newNode = node;
-            for(int i = 0; i<theNbSteps; i++) {
-              aTrsf2.Transforms( coord[0], coord[1], coord[2] );
-              newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
-              myLastCreatedNodes.Append(newNode);
-              listNewNodes.push_back( newNode );
-              srcNodes.Append( node );
-              aTrsf2.Transforms( coord[0], coord[1], coord[2] );
-              newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
-              myLastCreatedNodes.Append(newNode);
-              srcNodes.Append( node );
-              listNewNodes.push_back( newNode );
-            }
-          }
+        list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
+        if(listNewNodes.size()==theNbSteps) {
+        listNewNodes.clear();
+        // make new nodes
+        //gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
+        //double coord[3];
+        //aXYZ.Coord( coord[0], coord[1], coord[2] );
+        const SMDS_MeshNode * newNode = node;
+        if ( !isOnAxis ) {
+        for(int i = 0; i<theNbSteps; i++) {
+        aTrsf2.Transforms( coord[0], coord[1], coord[2] );
+        newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
+        cout<<"    3 AddNode:  "<<newNode;
+        myLastCreatedNodes.Append(newNode);
+        listNewNodes.push_back( newNode );
+        srcNodes.Append( node );
+        aTrsf2.Transforms( coord[0], coord[1], coord[2] );
+        newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
+        cout<<"    4 AddNode:  "<<newNode;
+        myLastCreatedNodes.Append(newNode);
+        srcNodes.Append( node );
+        listNewNodes.push_back( newNode );
         }
-      }
+        }
+        else {
+        listNewNodes.push_back( newNode );
+        }
+        }
+        }
+        }
+      */
       newNodesItVec.push_back( nIt );
     }
     // make new elements
@@ -3509,7 +3669,7 @@ SMESH_MeshEditor::RotationSweep(TIDSortedElemSet & theElems,
 
   if ( theMakeWalls )
     makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElems, theNbSteps, srcElems );
-  
+
   PGroupIDs newGroupIDs;
   if ( theMakeGroups )
     newGroupIDs = generateGroups( srcNodes, srcElems, "rotated");
@@ -3741,82 +3901,75 @@ SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet &  theElems,
   return newGroupIDs;
 }
 
-
+/*
 //=======================================================================
 //class    : SMESH_MeshEditor_PathPoint
 //purpose  : auxiliary class
 //=======================================================================
 class SMESH_MeshEditor_PathPoint {
 public:
-  SMESH_MeshEditor_PathPoint() {
-    myPnt.SetCoord(99., 99., 99.);
-    myTgt.SetCoord(1.,0.,0.);
-    myAngle=0.;
-    myPrm=0.;
-  }
-  void SetPnt(const gp_Pnt& aP3D){
-    myPnt=aP3D;
-  }
-  void SetTangent(const gp_Dir& aTgt){
-    myTgt=aTgt;
-  }
-  void SetAngle(const double& aBeta){
-    myAngle=aBeta;
-  }
-  void SetParameter(const double& aPrm){
-    myPrm=aPrm;
-  }
-  const gp_Pnt& Pnt()const{
-    return myPnt;
-  }
-  const gp_Dir& Tangent()const{
-    return myTgt;
-  }
-  double Angle()const{
-    return myAngle;
-  }
-  double Parameter()const{
-    return myPrm;
-  }
+SMESH_MeshEditor_PathPoint() {
+myPnt.SetCoord(99., 99., 99.);
+myTgt.SetCoord(1.,0.,0.);
+myAngle=0.;
+myPrm=0.;
+}
+void SetPnt(const gp_Pnt& aP3D){
+myPnt=aP3D;
+}
+void SetTangent(const gp_Dir& aTgt){
+myTgt=aTgt;
+}
+void SetAngle(const double& aBeta){
+myAngle=aBeta;
+}
+void SetParameter(const double& aPrm){
+myPrm=aPrm;
+}
+const gp_Pnt& Pnt()const{
+return myPnt;
+}
+const gp_Dir& Tangent()const{
+return myTgt;
+}
+double Angle()const{
+return myAngle;
+}
+double Parameter()const{
+return myPrm;
+}
 
 protected:
-  gp_Pnt myPnt;
-  gp_Dir myTgt;
-  double myAngle;
-  double myPrm;
+gp_Pnt myPnt;
+gp_Dir myTgt;
+double myAngle;
+double myPrm;
 };
+*/
 
 //=======================================================================
 //function : ExtrusionAlongTrack
 //purpose  :
 //=======================================================================
 SMESH_MeshEditor::Extrusion_Error
-  SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
-                                        SMESH_subMesh*       theTrack,
-                                        const SMDS_MeshNode* theN1,
-                                        const bool           theHasAngles,
-                                        list<double>&        theAngles,
-                                        const bool           theHasRefPoint,
-                                        const gp_Pnt&        theRefPoint,
-                                         const bool           theMakeGroups)
+SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
+                                       SMESH_subMesh*       theTrack,
+                                       const SMDS_MeshNode* theN1,
+                                       const bool           theHasAngles,
+                                       list<double>&        theAngles,
+                                       const bool           theLinearVariation,
+                                       const bool           theHasRefPoint,
+                                       const gp_Pnt&        theRefPoint,
+                                       const bool           theMakeGroups)
 {
   myLastCreatedElems.Clear();
   myLastCreatedNodes.Clear();
 
-  // source elements for each generated one
-  SMESH_SequenceOfElemPtr srcElems, srcNodes;
-
-  int j, aNbTP, aNbE, aNb;
-  double aT1, aT2, aT, aAngle, aX, aY, aZ;
+  int aNbE;
   std::list<double> aPrms;
-  std::list<double>::iterator aItD;
   TIDSortedElemSet::iterator itElem;
 
-  Standard_Real aTx1, aTx2, aL2, aTolVec, aTolVec2;
-  gp_Pnt aP3D, aV0;
-  gp_Vec aVec;
   gp_XYZ aGC;
-  Handle(Geom_Curve) aC3D;
   TopoDS_Edge aTrackEdge;
   TopoDS_Vertex aV1, aV2;
 
@@ -3825,11 +3978,6 @@ SMESH_MeshEditor::Extrusion_Error
   SMDSAbs_ElementType aTypeE;
 
   TNodeOfNodeListMap mapNewNodes;
-  TElemOfVecOfNnlmiMap mapElemNewNodes;
-  TElemOfElemListMap newElemsMap;
-
-  aTolVec=1.e-7;
-  aTolVec2=aTolVec*aTolVec;
 
   // 1. Check data
   aNbE = theElements.size();
@@ -3840,7 +3988,7 @@ SMESH_MeshEditor::Extrusion_Error
   // 1.1 Track Pattern
   ASSERT( theTrack );
 
-  SMESHDS_SubMesh* pSubMeshDS=theTrack->GetSubMeshDS();
+  SMESHDS_SubMesh* pSubMeshDS = theTrack->GetSubMeshDS();
 
   aItE = pSubMeshDS->GetElements();
   while ( aItE->more() ) {
@@ -3851,63 +3999,327 @@ SMESH_MeshEditor::Extrusion_Error
       return EXTR_PATH_NOT_EDGE;
   }
 
+  list<SMESH_MeshEditor_PathPoint> fullList;
+
   const TopoDS_Shape& aS = theTrack->GetSubShape();
-  // Sub shape for the Pattern must be an Edge
-  if ( aS.ShapeType() != TopAbs_EDGE )
+  // Sub shape for the Pattern must be an Edge or Wire
+  if( aS.ShapeType() == TopAbs_EDGE ) {
+    aTrackEdge = TopoDS::Edge( aS );
+    // the Edge must not be degenerated
+    if ( BRep_Tool::Degenerated( aTrackEdge ) )
+      return EXTR_BAD_PATH_SHAPE;
+    TopExp::Vertices( aTrackEdge, aV1, aV2 );
+    aItN = theTrack->GetFather()->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
+    const SMDS_MeshNode* aN1 = aItN->next();
+    aItN = theTrack->GetFather()->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
+    const SMDS_MeshNode* aN2 = aItN->next();
+    // starting node must be aN1 or aN2
+    if ( !( aN1 == theN1 || aN2 == theN1 ) )
+      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().get() );
+      double aT = pEPos->GetUParameter();
+      aPrms.push_back( aT );
+    }
+    //Extrusion_Error err =
+    MakeEdgePathPoints(aPrms, aTrackEdge, (aN1==theN1), fullList);
+  }
+  else if( aS.ShapeType() == TopAbs_WIRE ) {
+    list< SMESH_subMesh* > LSM;
+    TopTools_SequenceOfShape Edges;
+    SMESH_subMeshIteratorPtr itSM = theTrack->getDependsOnIterator(false,true);
+    while(itSM->more()) {
+      SMESH_subMesh* SM = itSM->next();
+      LSM.push_back(SM);
+      const TopoDS_Shape& aS = SM->GetSubShape();
+      Edges.Append(aS);
+    }
+    list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
+    int startNid = theN1->GetID();
+    TColStd_MapOfInteger UsedNums;
+    int NbEdges = Edges.Length();
+    int i = 1;
+    for(; i<=NbEdges; i++) {
+      int k = 0;
+      list< SMESH_subMesh* >::iterator itLSM = LSM.begin();
+      for(; itLSM!=LSM.end(); itLSM++) {
+        k++;
+        if(UsedNums.Contains(k)) continue;
+        aTrackEdge = TopoDS::Edge( Edges.Value(k) );
+        SMESH_subMesh* locTrack = *itLSM;
+        SMESHDS_SubMesh* locMeshDS = locTrack->GetSubMeshDS();
+        TopExp::Vertices( aTrackEdge, aV1, aV2 );
+        aItN = locTrack->GetFather()->GetSubMesh(aV1)->GetSubMeshDS()->GetNodes();
+        const SMDS_MeshNode* aN1 = aItN->next();
+        aItN = locTrack->GetFather()->GetSubMesh(aV2)->GetSubMeshDS()->GetNodes();
+        const SMDS_MeshNode* aN2 = aItN->next();
+        // starting node must be aN1 or aN2
+        if ( !( aN1->GetID() == startNid || aN2->GetID() == startNid ) ) continue;
+        // 2. Collect parameters on the track edge
+        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().get() );
+          double aT = pEPos->GetUParameter();
+          aPrms.push_back( aT );
+        }
+        list<SMESH_MeshEditor_PathPoint> LPP;
+        //Extrusion_Error err =
+        MakeEdgePathPoints(aPrms, aTrackEdge,(aN1->GetID()==startNid), LPP);
+        LLPPs.push_back(LPP);
+        UsedNums.Add(k);
+        // update startN for search following egde
+        if( aN1->GetID() == startNid ) startNid = aN2->GetID();
+        else startNid = aN1->GetID();
+        break;
+      }
+    }
+    list< list<SMESH_MeshEditor_PathPoint> >::iterator itLLPP = LLPPs.begin();
+    list<SMESH_MeshEditor_PathPoint> firstList = *itLLPP;
+    list<SMESH_MeshEditor_PathPoint>::iterator itPP = firstList.begin();
+    for(; itPP!=firstList.end(); itPP++) {
+      fullList.push_back( *itPP );
+    }
+    SMESH_MeshEditor_PathPoint PP1 = fullList.back();
+    fullList.pop_back();
+    itLLPP++;
+    for(; itLLPP!=LLPPs.end(); itLLPP++) {
+      list<SMESH_MeshEditor_PathPoint> currList = *itLLPP;
+      itPP = currList.begin();
+      SMESH_MeshEditor_PathPoint PP2 = currList.front();
+      gp_Dir D1 = PP1.Tangent();
+      gp_Dir D2 = PP2.Tangent();
+      gp_Dir Dnew( gp_Vec( (D1.X()+D2.X())/2, (D1.Y()+D2.Y())/2,
+                           (D1.Z()+D2.Z())/2 ) );
+      PP1.SetTangent(Dnew);
+      fullList.push_back(PP1);
+      itPP++;
+      for(; itPP!=firstList.end(); itPP++) {
+        fullList.push_back( *itPP );
+      }
+      PP1 = fullList.back();
+      fullList.pop_back();
+    }
+    // if wire not closed
+    fullList.push_back(PP1);
+    // else ???
+  }
+  else {
     return EXTR_BAD_PATH_SHAPE;
+  }
 
-  aTrackEdge = TopoDS::Edge( aS );
-  // the Edge must not be degenerated
-  if ( BRep_Tool::Degenerated( aTrackEdge ) )
-    return EXTR_BAD_PATH_SHAPE;
+  return MakeExtrElements(theElements, fullList, theHasAngles, theAngles, theLinearVariation,
+                          theHasRefPoint, theRefPoint, theMakeGroups);
+}
 
-  TopExp::Vertices( aTrackEdge, aV1, aV2 );
-  aT1=BRep_Tool::Parameter( aV1, aTrackEdge );
-  aT2=BRep_Tool::Parameter( aV2, aTrackEdge );
 
-  aItN = theTrack->GetFather()->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
-  const SMDS_MeshNode* aN1 = aItN->next();
+//=======================================================================
+//function : ExtrusionAlongTrack
+//purpose  :
+//=======================================================================
+SMESH_MeshEditor::Extrusion_Error
+SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
+                                       SMESH_Mesh*          theTrack,
+                                       const SMDS_MeshNode* theN1,
+                                       const bool           theHasAngles,
+                                       list<double>&        theAngles,
+                                       const bool           theLinearVariation,
+                                       const bool           theHasRefPoint,
+                                       const gp_Pnt&        theRefPoint,
+                                       const bool           theMakeGroups)
+{
+  myLastCreatedElems.Clear();
+  myLastCreatedNodes.Clear();
+
+  int aNbE;
+  std::list<double> aPrms;
+  TIDSortedElemSet::iterator itElem;
 
-  aItN = theTrack->GetFather()->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
-  const SMDS_MeshNode* aN2 = aItN->next();
+  gp_XYZ aGC;
+  TopoDS_Edge aTrackEdge;
+  TopoDS_Vertex aV1, aV2;
 
-  // starting node must be aN1 or aN2
-  if ( !( aN1 == theN1 || aN2 == theN1 ) )
-    return EXTR_BAD_STARTING_NODE;
+  SMDS_ElemIteratorPtr aItE;
+  SMDS_NodeIteratorPtr aItN;
+  SMDSAbs_ElementType aTypeE;
 
-  aNbTP = pSubMeshDS->NbNodes() + 2;
+  TNodeOfNodeListMap mapNewNodes;
 
-  // 1.2. Angles
-  vector<double> aAngles( aNbTP );
+  // 1. Check data
+  aNbE = theElements.size();
+  // nothing to do
+  if ( !aNbE )
+    return EXTR_NO_ELEMENTS;
 
-  for ( j=0; j < aNbTP; ++j ) {
-    aAngles[j] = 0.;
+  // 1.1 Track Pattern
+  ASSERT( theTrack );
+
+  SMESHDS_Mesh* pMeshDS = theTrack->GetMeshDS();
+
+  aItE = pMeshDS->elementsIterator();
+  while ( aItE->more() ) {
+    const SMDS_MeshElement* pE = aItE->next();
+    aTypeE = pE->GetType();
+    // Pattern must contain links only
+    if ( aTypeE != SMDSAbs_Edge )
+      return EXTR_PATH_NOT_EDGE;
   }
 
-  if ( theHasAngles ) {
-    aItD = theAngles.begin();
-    for ( j=1; (aItD != theAngles.end()) && (j<aNbTP); ++aItD, ++j ) {
-      aAngle = *aItD;
-      aAngles[j] = aAngle;
+  list<SMESH_MeshEditor_PathPoint> fullList;
+
+  const TopoDS_Shape& aS = theTrack->GetShapeToMesh();
+  // Sub shape for the Pattern must be an Edge or Wire
+  if( aS.ShapeType() == TopAbs_EDGE ) {
+    aTrackEdge = TopoDS::Edge( aS );
+    // the Edge must not be degenerated
+    if ( BRep_Tool::Degenerated( aTrackEdge ) )
+      return EXTR_BAD_PATH_SHAPE;
+    TopExp::Vertices( aTrackEdge, aV1, aV2 );
+    aItN = theTrack->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
+    const SMDS_MeshNode* aN1 = aItN->next();
+    aItN = theTrack->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
+    const SMDS_MeshNode* aN2 = aItN->next();
+    // starting node must be aN1 or aN2
+    if ( !( aN1 == theN1 || aN2 == theN1 ) )
+      return EXTR_BAD_STARTING_NODE;
+    aItN = pMeshDS->nodesIterator();
+    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().get() );
+      double aT = pEPos->GetUParameter();
+      aPrms.push_back( aT );
+    }
+    //Extrusion_Error err =
+    MakeEdgePathPoints(aPrms, aTrackEdge, (aN1==theN1), fullList);
+  }
+  else if( aS.ShapeType() == TopAbs_WIRE ) {
+    list< SMESH_subMesh* > LSM;
+    TopTools_SequenceOfShape Edges;
+    TopExp_Explorer eExp(aS, TopAbs_EDGE);
+    for(; eExp.More(); eExp.Next()) {
+      TopoDS_Edge E = TopoDS::Edge( eExp.Current() );
+      if( BRep_Tool::Degenerated(E) ) continue;
+      SMESH_subMesh* SM = theTrack->GetSubMesh(E);
+      if(SM) {
+        LSM.push_back(SM);
+        Edges.Append(E);
+      }
+    }
+    list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
+    int startNid = theN1->GetID();
+    TColStd_MapOfInteger UsedNums;
+    int NbEdges = Edges.Length();
+    int i = 1;
+    for(; i<=NbEdges; i++) {
+      int k = 0;
+      list< SMESH_subMesh* >::iterator itLSM = LSM.begin();
+      for(; itLSM!=LSM.end(); itLSM++) {
+        k++;
+        if(UsedNums.Contains(k)) continue;
+        aTrackEdge = TopoDS::Edge( Edges.Value(k) );
+        SMESH_subMesh* locTrack = *itLSM;
+        SMESHDS_SubMesh* locMeshDS = locTrack->GetSubMeshDS();
+        TopExp::Vertices( aTrackEdge, aV1, aV2 );
+        aItN = locTrack->GetFather()->GetSubMesh(aV1)->GetSubMeshDS()->GetNodes();
+        const SMDS_MeshNode* aN1 = aItN->next();
+        aItN = locTrack->GetFather()->GetSubMesh(aV2)->GetSubMeshDS()->GetNodes();
+        const SMDS_MeshNode* aN2 = aItN->next();
+        // starting node must be aN1 or aN2
+        if ( !( aN1->GetID() == startNid || aN2->GetID() == startNid ) ) continue;
+        // 2. Collect parameters on the track edge
+        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().get() );
+          double aT = pEPos->GetUParameter();
+          aPrms.push_back( aT );
+        }
+        list<SMESH_MeshEditor_PathPoint> LPP;
+        //Extrusion_Error err =
+        MakeEdgePathPoints(aPrms, aTrackEdge,(aN1->GetID()==startNid), LPP);
+        LLPPs.push_back(LPP);
+        UsedNums.Add(k);
+        // update startN for search following egde
+        if( aN1->GetID() == startNid ) startNid = aN2->GetID();
+        else startNid = aN1->GetID();
+        break;
+      }
     }
+    list< list<SMESH_MeshEditor_PathPoint> >::iterator itLLPP = LLPPs.begin();
+    list<SMESH_MeshEditor_PathPoint> firstList = *itLLPP;
+    list<SMESH_MeshEditor_PathPoint>::iterator itPP = firstList.begin();
+    for(; itPP!=firstList.end(); itPP++) {
+      fullList.push_back( *itPP );
+    }
+    SMESH_MeshEditor_PathPoint PP1 = fullList.back();
+    fullList.pop_back();
+    itLLPP++;
+    for(; itLLPP!=LLPPs.end(); itLLPP++) {
+      list<SMESH_MeshEditor_PathPoint> currList = *itLLPP;
+      itPP = currList.begin();
+      SMESH_MeshEditor_PathPoint PP2 = currList.front();
+      gp_Pnt P1 = PP1.Pnt();
+      //cout<<"    PP1: Pnt("<<P1.X()<<","<<P1.Y()<<","<<P1.Z()<<")"<<endl;
+      gp_Pnt P2 = PP2.Pnt();
+      gp_Dir D1 = PP1.Tangent();
+      gp_Dir D2 = PP2.Tangent();
+      gp_Dir Dnew( gp_Vec( (D1.X()+D2.X())/2, (D1.Y()+D2.Y())/2,
+                           (D1.Z()+D2.Z())/2 ) );
+      PP1.SetTangent(Dnew);
+      fullList.push_back(PP1);
+      itPP++;
+      for(; itPP!=currList.end(); itPP++) {
+        fullList.push_back( *itPP );
+      }
+      PP1 = fullList.back();
+      fullList.pop_back();
+    }
+    // if wire not closed
+    fullList.push_back(PP1);
+    // else ???
+  }
+  else {
+    return EXTR_BAD_PATH_SHAPE;
   }
 
-  // 2. Collect parameters on the track edge
-  aPrms.push_back( aT1 );
-  aPrms.push_back( aT2 );
+  return MakeExtrElements(theElements, fullList, theHasAngles, theAngles, theLinearVariation,
+                          theHasRefPoint, theRefPoint, theMakeGroups);
+}
 
-  aItN = pSubMeshDS->GetNodes();
-  while ( aItN->more() ) {
-    const SMDS_MeshNode* pNode = aItN->next();
-    const SMDS_EdgePosition* pEPos =
-      static_cast<const SMDS_EdgePosition*>( pNode->GetPosition().get() );
-    aT = pEPos->GetUParameter();
-    aPrms.push_back( aT );
-  }
 
+//=======================================================================
+//function : MakeEdgePathPoints
+//purpose  : auxilary for ExtrusionAlongTrack
+//=======================================================================
+SMESH_MeshEditor::Extrusion_Error
+SMESH_MeshEditor::MakeEdgePathPoints(std::list<double>& aPrms,
+                                     const TopoDS_Edge& aTrackEdge,
+                                     bool FirstIsStart,
+                                     list<SMESH_MeshEditor_PathPoint>& LPP)
+{
+  Standard_Real aTx1, aTx2, aL2, aTolVec, aTolVec2;
+  aTolVec=1.e-7;
+  aTolVec2=aTolVec*aTolVec;
+  double aT1, aT2;
+  TopoDS_Vertex aV1, aV2;
+  TopExp::Vertices( aTrackEdge, aV1, aV2 );
+  aT1=BRep_Tool::Parameter( aV1, aTrackEdge );
+  aT2=BRep_Tool::Parameter( aV2, aTrackEdge );
+  // 2. Collect parameters on the track edge
+  aPrms.push_front( aT1 );
+  aPrms.push_back( aT2 );
   // sort parameters
   aPrms.sort();
-  if ( aN1 == theN1 ) {
+  if( FirstIsStart ) {
     if ( aT1 > aT2 ) {
       aPrms.reverse();
     }
@@ -3917,33 +4329,86 @@ SMESH_MeshEditor::Extrusion_Error
       aPrms.reverse();
     }
   }
-
   // 3. Path Points
   SMESH_MeshEditor_PathPoint aPP;
-  vector<SMESH_MeshEditor_PathPoint> aPPs( aNbTP );
-  //
-  aC3D = BRep_Tool::Curve( aTrackEdge, aTx1, aTx2 );
-  //
-  aItD = aPrms.begin();
-  for ( j=0; aItD != aPrms.end(); ++aItD, ++j ) {
-    aT = *aItD;
+  Handle(Geom_Curve) aC3D = BRep_Tool::Curve( aTrackEdge, aTx1, aTx2 );
+  std::list<double>::iterator aItD = aPrms.begin();
+  for(; aItD != aPrms.end(); ++aItD) {
+    double aT = *aItD;
+    gp_Pnt aP3D;
+    gp_Vec aVec;
     aC3D->D1( aT, aP3D, aVec );
     aL2 = aVec.SquareMagnitude();
     if ( aL2 < aTolVec2 )
       return EXTR_CANT_GET_TANGENT;
-
     gp_Dir aTgt( aVec );
-    aAngle = aAngles[j];
-
     aPP.SetPnt( aP3D );
     aPP.SetTangent( aTgt );
-    aPP.SetAngle( aAngle );
     aPP.SetParameter( aT );
-    aPPs[j]=aPP;
+    LPP.push_back(aPP);
+  }
+  return EXTR_OK;
+}
+
+
+//=======================================================================
+//function : MakeExtrElements
+//purpose  : auxilary for ExtrusionAlongTrack
+//=======================================================================
+SMESH_MeshEditor::Extrusion_Error
+SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet&  theElements,
+                                   list<SMESH_MeshEditor_PathPoint>& fullList,
+                                   const bool theHasAngles,
+                                   list<double>& theAngles,
+                                   const bool theLinearVariation,
+                                   const bool theHasRefPoint,
+                                   const gp_Pnt& theRefPoint,
+                                   const bool theMakeGroups)
+{
+  //cout<<"MakeExtrElements  fullList.size() = "<<fullList.size()<<endl;
+  int aNbTP = fullList.size();
+  vector<SMESH_MeshEditor_PathPoint> aPPs(aNbTP);
+  // Angles
+  if( theHasAngles && theAngles.size()>0 && theLinearVariation ) {
+    LinearAngleVariation(aNbTP-1, theAngles);
+  }
+  vector<double> aAngles( aNbTP );
+  int j = 0;
+  for(; j<aNbTP; ++j) {
+    aAngles[j] = 0.;
+  }
+  if ( theHasAngles ) {
+    double anAngle;;
+    std::list<double>::iterator aItD = theAngles.begin();
+    for ( j=1; (aItD != theAngles.end()) && (j<aNbTP); ++aItD, ++j ) {
+      anAngle = *aItD;
+      aAngles[j] = anAngle;
+    }
+  }
+  // fill vector of path points with angles
+  //aPPs.resize(fullList.size());
+  j = -1;
+  list<SMESH_MeshEditor_PathPoint>::iterator itPP = fullList.begin();
+  for(; itPP!=fullList.end(); itPP++) {
+    j++;
+    SMESH_MeshEditor_PathPoint PP = *itPP;
+    PP.SetAngle(aAngles[j]);
+    aPPs[j] = PP;
   }
 
+  TNodeOfNodeListMap mapNewNodes;
+  TElemOfVecOfNnlmiMap mapElemNewNodes;
+  TElemOfElemListMap newElemsMap;
+  TIDSortedElemSet::iterator itElem;
+  double aX, aY, aZ;
+  int aNb;
+  SMDSAbs_ElementType aTypeE;
+  // source elements for each generated one
+  SMESH_SequenceOfElemPtr srcElems, srcNodes;
+
   // 3. Center of rotation aV0
-  aV0 = theRefPoint;
+  gp_Pnt aV0 = theRefPoint;
+  gp_XYZ aGC;
   if ( !theHasRefPoint ) {
     aNb = 0;
     aGC.SetCoord( 0.,0.,0. );
@@ -3954,19 +4419,19 @@ SMESH_MeshEditor::Extrusion_Error
 
       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
       while ( itN->more() ) {
-       const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( itN->next() );
-       aX = node->X();
-       aY = node->Y();
-       aZ = node->Z();
-
-       if ( mapNewNodes.find( node ) == mapNewNodes.end() ) {
-         list<const SMDS_MeshNode*> aLNx;
-         mapNewNodes[node] = aLNx;
-         //
-         gp_XYZ aXYZ( aX, aY, aZ );
-         aGC += aXYZ;
-         ++aNb;
-       }
+        const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( itN->next() );
+        aX = node->X();
+        aY = node->Y();
+        aZ = node->Z();
+
+        if ( mapNewNodes.find( node ) == mapNewNodes.end() ) {
+          list<const SMDS_MeshNode*> aLNx;
+          mapNewNodes[node] = aLNx;
+          //
+          gp_XYZ aXYZ( aX, aY, aZ );
+          aGC += aXYZ;
+          ++aNb;
+        }
       }
     }
     aGC /= aNb;
@@ -3995,65 +4460,66 @@ SMESH_MeshEditor::Extrusion_Error
       ++nodeIndex;
       // check if a node has been already processed
       const SMDS_MeshNode* node =
-       static_cast<const SMDS_MeshNode*>( itN->next() );
+        static_cast<const SMDS_MeshNode*>( itN->next() );
       TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
       if ( nIt == mapNewNodes.end() ) {
         nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
         list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
 
-       // make new nodes
-       aX = node->X();  aY = node->Y(); aZ = node->Z();
-
-       Standard_Real aAngle1x, aAngleT1T0, aTolAng;
-       gp_Pnt aP0x, aP1x, aPN0, aPN1, aV0x, aV1x;
-       gp_Ax1 anAx1, anAxT1T0;
-       gp_Dir aDT1x, aDT0x, aDT1T0;
-
-       aTolAng=1.e-4;
-
-       aV0x = aV0;
-       aPN0.SetCoord(aX, aY, aZ);
+        // make new nodes
+        aX = node->X();  aY = node->Y(); aZ = node->Z();
+
+        Standard_Real aAngle1x, aAngleT1T0, aTolAng;
+        gp_Pnt aP0x, aP1x, aPN0, aPN1, aV0x, aV1x;
+        gp_Ax1 anAx1, anAxT1T0;
+        gp_Dir aDT1x, aDT0x, aDT1T0;
+
+        aTolAng=1.e-4;
+
+        aV0x = aV0;
+        aPN0.SetCoord(aX, aY, aZ);
+
+        const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0];
+        aP0x = aPP0.Pnt();
+        aDT0x= aPP0.Tangent();
+        //cout<<"j = 0   PP: Pnt("<<aP0x.X()<<","<<aP0x.Y()<<","<<aP0x.Z()<<")"<<endl;
+
+        for ( j = 1; j < aNbTP; ++j ) {
+          const SMESH_MeshEditor_PathPoint& aPP1 = aPPs[j];
+          aP1x = aPP1.Pnt();
+          aDT1x = aPP1.Tangent();
+          aAngle1x = aPP1.Angle();
+
+          gp_Trsf aTrsf, aTrsfRot, aTrsfRotT1T0;
+          // Translation
+          gp_Vec aV01x( aP0x, aP1x );
+          aTrsf.SetTranslation( aV01x );
+
+          // traslated point
+          aV1x = aV0x.Transformed( aTrsf );
+          aPN1 = aPN0.Transformed( aTrsf );
+
+          // rotation 1 [ T1,T0 ]
+          aAngleT1T0=-aDT1x.Angle( aDT0x );
+          if (fabs(aAngleT1T0) > aTolAng) {
+            aDT1T0=aDT1x^aDT0x;
+            anAxT1T0.SetLocation( aV1x );
+            anAxT1T0.SetDirection( aDT1T0 );
+            aTrsfRotT1T0.SetRotation( anAxT1T0, aAngleT1T0 );
+
+            aPN1 = aPN1.Transformed( aTrsfRotT1T0 );
+          }
 
-       const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0];
-       aP0x = aPP0.Pnt();
-       aDT0x= aPP0.Tangent();
+          // rotation 2
+          if ( theHasAngles ) {
+            anAx1.SetLocation( aV1x );
+            anAx1.SetDirection( aDT1x );
+            aTrsfRot.SetRotation( anAx1, aAngle1x );
 
-       for ( j = 1; j < aNbTP; ++j ) {
-         const SMESH_MeshEditor_PathPoint& aPP1 = aPPs[j];
-         aP1x = aPP1.Pnt();
-         aDT1x = aPP1.Tangent();
-         aAngle1x = aPP1.Angle();
+            aPN1 = aPN1.Transformed( aTrsfRot );
+          }
 
-         gp_Trsf aTrsf, aTrsfRot, aTrsfRotT1T0;
-         // Translation
-         gp_Vec aV01x( aP0x, aP1x );
-         aTrsf.SetTranslation( aV01x );
-
-         // traslated point
-         aV1x = aV0x.Transformed( aTrsf );
-         aPN1 = aPN0.Transformed( aTrsf );
-
-         // rotation 1 [ T1,T0 ]
-         aAngleT1T0=-aDT1x.Angle( aDT0x );
-         if (fabs(aAngleT1T0) > aTolAng) {
-           aDT1T0=aDT1x^aDT0x;
-           anAxT1T0.SetLocation( aV1x );
-           anAxT1T0.SetDirection( aDT1T0 );
-           aTrsfRotT1T0.SetRotation( anAxT1T0, aAngleT1T0 );
-
-           aPN1 = aPN1.Transformed( aTrsfRotT1T0 );
-         }
-
-         // rotation 2
-         if ( theHasAngles ) {
-           anAx1.SetLocation( aV1x );
-           anAx1.SetDirection( aDT1x );
-           aTrsfRot.SetRotation( anAx1, aAngle1x );
-
-           aPN1 = aPN1.Transformed( aTrsfRot );
-         }
-
-         // make new node
+          // make new node
           if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
             // create additional node
             double x = ( aPN1.X() + aPN0.X() )/2.;
@@ -4064,19 +4530,19 @@ SMESH_MeshEditor::Extrusion_Error
             srcNodes.Append( node );
             listNewNodes.push_back( newNode );
           }
-         aX = aPN1.X();
-         aY = aPN1.Y();
-         aZ = aPN1.Z();
-         const SMDS_MeshNode* newNode = aMesh->AddNode( aX, aY, aZ );
+          aX = aPN1.X();
+          aY = aPN1.Y();
+          aZ = aPN1.Z();
+          const SMDS_MeshNode* newNode = aMesh->AddNode( aX, aY, aZ );
           myLastCreatedNodes.Append(newNode);
           srcNodes.Append( node );
-         listNewNodes.push_back( newNode );
+          listNewNodes.push_back( newNode );
 
-         aPN0 = aPN1;
-         aP0x = aP1x;
-         aV0x = aV1x;
-         aDT0x = aDT1x;
-       }
+          aPN0 = aPN1;
+          aP0x = aP1x;
+          aV0x = aV1x;
+          aDT0x = aDT1x;
+        }
       }
 
       else {
@@ -4125,6 +4591,56 @@ SMESH_MeshEditor::Extrusion_Error
   return EXTR_OK;
 }
 
+
+//=======================================================================
+//function : LinearAngleVariation
+//purpose  : auxilary for ExtrusionAlongTrack
+//=======================================================================
+void SMESH_MeshEditor::LinearAngleVariation(const int nbSteps,
+                                            list<double>& Angles)
+{
+  int nbAngles = Angles.size();
+  if( nbSteps > nbAngles ) {
+    vector<double> theAngles(nbAngles);
+    list<double>::iterator it = Angles.begin();
+    int i = -1;
+    for(; it!=Angles.end(); it++) {
+      i++;
+      theAngles[i] = (*it);
+    }
+    list<double> res;
+    double rAn2St = double( nbAngles ) / double( nbSteps );
+    double angPrev = 0, angle;
+    for ( int iSt = 0; iSt < nbSteps; ++iSt ) {
+      double angCur = rAn2St * ( iSt+1 );
+      double angCurFloor  = floor( angCur );
+      double angPrevFloor = floor( angPrev );
+      if ( angPrevFloor == angCurFloor )
+        angle = rAn2St * theAngles[ int( angCurFloor ) ];
+      else {
+        int iP = int( angPrevFloor );
+        double angPrevCeil = ceil(angPrev);
+        angle = ( angPrevCeil - angPrev ) * theAngles[ iP ];
+
+        int iC = int( angCurFloor );
+        if ( iC < nbAngles )
+          angle += ( angCur - angCurFloor ) * theAngles[ iC ];
+
+        iP = int( angPrevCeil );
+        while ( iC-- > iP )
+          angle += theAngles[ iC ];
+      }
+      res.push_back(angle);
+      angPrev = angCur;
+    }
+    Angles.clear();
+    it = res.begin();
+    for(; it!=res.end(); it++)
+      Angles.push_back( *it );
+  }
+}
+
+
 //=======================================================================
 //function : Transform
 //purpose  :
@@ -4166,7 +4682,7 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
   SMESH_MeshEditor targetMeshEditor( theTargetMesh );
   SMESHDS_Mesh* aTgtMesh = theTargetMesh ? theTargetMesh->GetMeshDS() : 0;
   SMESHDS_Mesh* aMesh    = GetMeshDS();
-  
+
 
   // map old node to new one
   TNodeNodeMap nodeMap;
@@ -4248,7 +4764,7 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
     REV_FACE    = 3,
     REV_HEXA    = 4,  //  = nbNodes - 4
     FORWARD     = 5
-    };
+  };
   int index[][8] = {
     { 2, 1, 0, 3, 4, 0, 0, 0 },  // REV_TETRA
     { 2, 1, 0, 3, 4, 0, 0, 0 },  // REV_PYRAMID
@@ -4349,18 +4865,18 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
           }
         }
         break;
-    default:;
+      default:;
+      }
+      continue;
     }
-    continue;
-  }
 
-  // Regular elements
-  int* i = index[ FORWARD ];
-  if ( needReverse && nbNodes > 2) // reverse mirrored faces and volumes
-    if ( elemType == SMDSAbs_Face )
-      i = index[ REV_FACE ];
-    else
-      i = index[ nbNodes - 4 ];
+    // Regular elements
+    int* i = index[ FORWARD ];
+    if ( needReverse && nbNodes > 2) // reverse mirrored faces and volumes
+      if ( elemType == SMDSAbs_Face )
+        i = index[ REV_FACE ];
+      else
+        i = index[ nbNodes - 4 ];
 
     if(elem->IsQuadratic()) {
       static int anIds[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19};
@@ -4567,12 +5083,13 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
   return newGroupIDs;
 }
 
-//=======================================================================
-//function : FindCoincidentNodes
-//purpose  : Return list of group of nodes close to each other within theTolerance
-//           Search among theNodes or in the whole mesh if theNodes is empty using
-//           an Octree algorithm
-//=======================================================================
+//================================================================================
+/*!
+ * \brief Return list of group of nodes close to each other within theTolerance
+ *        Search among theNodes or in the whole mesh if theNodes is empty using
+ *        an Octree algorithm
+ */
+//================================================================================
 
 void SMESH_MeshEditor::FindCoincidentNodes (set<const SMDS_MeshNode*> & theNodes,
                                             const double                theTolerance,
@@ -4590,10 +5107,11 @@ void SMESH_MeshEditor::FindCoincidentNodes (set<const SMDS_MeshNode*> & theNodes
   }
   else
     nodes=theNodes;
-  SMESH_OctreeNode::FindCoincidentNodes ( nodes, &theGroupsOfNodes, theTolerance);
 
+  SMESH_OctreeNode::FindCoincidentNodes ( nodes, &theGroupsOfNodes, theTolerance);
 }
 
+
 //=======================================================================
 /*!
  * \brief Implementation of search for the node closest to point
@@ -4602,11 +5120,14 @@ void SMESH_MeshEditor::FindCoincidentNodes (set<const SMDS_MeshNode*> & theNodes
 
 struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
 {
+  //---------------------------------------------------------------------
   /*!
    * \brief Constructor
    */
   SMESH_NodeSearcherImpl( const SMESHDS_Mesh* theMesh )
   {
+    myMesh = ( SMESHDS_Mesh* ) theMesh;
+
     set<const SMDS_MeshNode*> nodes;
     if ( theMesh ) {
       SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator();
@@ -4614,19 +5135,43 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
         nodes.insert( nodes.end(), nIt->next() );
     }
     myOctreeNode = new SMESH_OctreeNode(nodes) ;
+
+    // get max size of a leaf box
+    SMESH_OctreeNode* tree = myOctreeNode;
+    while ( !tree->isLeaf() )
+    {
+      SMESH_OctreeNodeIteratorPtr cIt = tree->GetChildrenIterator();
+      if ( cIt->more() )
+        tree = cIt->next();
+    }
+    myHalfLeafSize = tree->maxSize() / 2.;
   }
+
+  //---------------------------------------------------------------------
+  /*!
+   * \brief Move node and update myOctreeNode accordingly
+   */
+  void MoveNode( const SMDS_MeshNode* node, const gp_Pnt& toPnt )
+  {
+    myOctreeNode->UpdateByMoveNode( node, toPnt );
+    myMesh->MoveNode( node, toPnt.X(), toPnt.Y(), toPnt.Z() );
+  }
+
+  //---------------------------------------------------------------------
   /*!
    * \brief Do it's job
    */
   const SMDS_MeshNode* FindClosestTo( const gp_Pnt& thePnt )
   {
     SMDS_MeshNode tgtNode( thePnt.X(), thePnt.Y(), thePnt.Z() );
+    map<double, const SMDS_MeshNode*> dist2Nodes;
+    myOctreeNode->NodesAround( &tgtNode, dist2Nodes, myHalfLeafSize );
+    if ( !dist2Nodes.empty() )
+      return dist2Nodes.begin()->second;
     list<const SMDS_MeshNode*> nodes;
-    const double precision = 1e-6;
-    myOctreeNode->NodesAround( &tgtNode, &nodes, precision );
+    //myOctreeNode->NodesAround( &tgtNode, &nodes, myHalfLeafSize );
 
     double minSqDist = DBL_MAX;
-    Bnd_B3d box;
     if ( nodes.empty() )  // get all nodes of OctreeNode's closest to thePnt
     {
       // sort leafs by their distance from thePnt
@@ -4635,20 +5180,25 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
       list< SMESH_OctreeNode* > treeList;
       list< SMESH_OctreeNode* >::iterator trIt;
       treeList.push_back( myOctreeNode );
+
+      SMDS_MeshNode pointNode( thePnt.X(), thePnt.Y(), thePnt.Z() );
       for ( trIt = treeList.begin(); trIt != treeList.end(); ++trIt)
       {
         SMESH_OctreeNode* tree = *trIt;
-        if ( !tree->isLeaf() ) { // put children to the queue
+        if ( !tree->isLeaf() ) // put children to the queue
+        {
+          if ( !tree->isInside( &pointNode, myHalfLeafSize )) continue;
           SMESH_OctreeNodeIteratorPtr cIt = tree->GetChildrenIterator();
           while ( cIt->more() )
             treeList.push_back( cIt->next() );
         }
-        else if ( tree->NbNodes() ) { // put tree to treeMap
-          tree->getBox( box );
+        else if ( tree->NbNodes() ) // put a tree to the treeMap
+        {
+          const Bnd_B3d& box = tree->getBox();
           double sqDist = thePnt.SquareDistance( 0.5 * ( box.CornerMin() + box.CornerMax() ));
           pair<TDistTreeMap::iterator,bool> it_in = treeMap.insert( make_pair( sqDist, tree ));
           if ( !it_in.second ) // not unique distance to box center
-            treeMap.insert( it_in.first, make_pair( sqDist - 1e-13*treeMap.size(), tree ));
+            treeMap.insert( it_in.first, make_pair( sqDist + 1e-13*treeMap.size(), tree ));
         }
       }
       // find distance after which there is no sense to check tree's
@@ -4656,7 +5206,7 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
       TDistTreeMap::iterator sqDist_tree = treeMap.begin();
       if ( treeMap.size() > 5 ) {
         SMESH_OctreeNode* closestTree = sqDist_tree->second;
-        closestTree->getBox( box );
+        const Bnd_B3d& box = closestTree->getBox();
         double limit = sqrt( sqDist_tree->first ) + sqrt ( box.SquareExtent() );
         sqLimit = limit * limit;
       }
@@ -4681,12 +5231,23 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
     }
     return closestNode;
   }
+
+  //---------------------------------------------------------------------
   /*!
    * \brief Destructor
    */
   ~SMESH_NodeSearcherImpl() { delete myOctreeNode; }
+
+  //---------------------------------------------------------------------
+  /*!
+   * \brief Return the node tree
+   */
+  const SMESH_OctreeNode* getTree() const { return myOctreeNode; }
+
 private:
   SMESH_OctreeNode* myOctreeNode;
+  SMESHDS_Mesh*     myMesh;
+  double            myHalfLeafSize; // max size of a leaf box
 };
 
 //=======================================================================
@@ -4700,6 +5261,404 @@ SMESH_NodeSearcher* SMESH_MeshEditor::GetNodeSearcher()
   return new SMESH_NodeSearcherImpl( GetMeshDS() );
 }
 
+// ========================================================================
+namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
+{
+  const int MaxNbElemsInLeaf = 10; // maximal number of elements in a leaf of tree
+  const int MaxLevel         = 7;  // maximal tree height -> nb terminal boxes: 8^7 = 2097152
+  const double NodeRadius = 1e-9;  // to enlarge bnd box of element
+
+  //=======================================================================
+  /*!
+   * \brief Octal tree of bounding boxes of elements
+   */
+  //=======================================================================
+
+  class ElementBndBoxTree : public SMESH_Octree
+  {
+  public:
+
+    ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType);
+    void getElementsNearPoint( const gp_Pnt& point, TIDSortedElemSet& foundElems);
+    ~ElementBndBoxTree();
+
+  protected:
+    ElementBndBoxTree() {}
+    SMESH_Octree* allocateOctreeChild() const { return new ElementBndBoxTree; }
+    void buildChildrenData();
+    Bnd_B3d* buildRootBox();
+  private:
+    //!< Bounding box of element
+    struct ElementBox : public Bnd_B3d
+    {
+      const SMDS_MeshElement* _element;
+      int                     _refCount; // an ElementBox can be included in several tree branches
+      ElementBox(const SMDS_MeshElement* elem);
+    };
+    vector< ElementBox* > _elements;
+  };
+
+  //================================================================================
+  /*!
+   * \brief ElementBndBoxTree creation
+   */
+  //================================================================================
+
+  ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType)
+    :SMESH_Octree( new SMESH_Octree::Limit( MaxLevel, /*minSize=*/0. ))
+  {
+    int nbElems = mesh.GetMeshInfo().NbElements( elemType );
+    _elements.reserve( nbElems );
+
+    SMDS_ElemIteratorPtr elemIt = mesh.elementsIterator( elemType );
+    while ( elemIt->more() )
+      _elements.push_back( new ElementBox( elemIt->next() ));
+
+    if ( _elements.size() > MaxNbElemsInLeaf )
+      compute();
+    else
+      myIsLeaf = true;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Destructor
+   */
+  //================================================================================
+
+  ElementBndBoxTree::~ElementBndBoxTree()
+  {
+    for ( int i = 0; i < _elements.size(); ++i )
+      if ( --_elements[i]->_refCount <= 0 )
+        delete _elements[i];
+  }
+
+  //================================================================================
+  /*!
+   * \brief Return the maximal box
+   */
+  //================================================================================
+
+  Bnd_B3d* ElementBndBoxTree::buildRootBox()
+  {
+    Bnd_B3d* box = new Bnd_B3d;
+    for ( int i = 0; i < _elements.size(); ++i )
+      box->Add( *_elements[i] );
+    return box;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Redistrubute element boxes among children
+   */
+  //================================================================================
+
+  void ElementBndBoxTree::buildChildrenData()
+  {
+    for ( int i = 0; i < _elements.size(); ++i )
+    {
+      for (int j = 0; j < 8; j++)
+      {
+        if ( !_elements[i]->IsOut( myChildren[j]->getBox() ))
+        {
+          _elements[i]->_refCount++;
+          ((ElementBndBoxTree*)myChildren[j])->_elements.push_back( _elements[i]);
+        }
+      }
+      _elements[i]->_refCount--;
+    }
+    _elements.clear();
+
+    for (int j = 0; j < 8; j++)
+    {
+      ElementBndBoxTree* child = static_cast<ElementBndBoxTree*>( myChildren[j]);
+      if ( child->_elements.size() <= MaxNbElemsInLeaf )
+        child->myIsLeaf = true;
+
+      if ( child->_elements.capacity() - child->_elements.size() > 1000 )
+        child->_elements.resize( child->_elements.size() ); // compact
+    }
+  }
+
+  //================================================================================
+  /*!
+   * \brief Return elements which can include the point
+   */
+  //================================================================================
+
+  void ElementBndBoxTree::getElementsNearPoint( const gp_Pnt&     point,
+                                                TIDSortedElemSet& foundElems)
+  {
+    if ( level() && getBox().IsOut( point.XYZ() ))
+      return;
+
+    if ( isLeaf() )
+    {
+      for ( int i = 0; i < _elements.size(); ++i )
+        if ( !_elements[i]->IsOut( point.XYZ() ))
+          foundElems.insert( _elements[i]->_element );
+    }
+    else
+    {
+      for (int i = 0; i < 8; i++)
+        ((ElementBndBoxTree*) myChildren[i])->getElementsNearPoint( point, foundElems );
+    }
+  }
+
+  //================================================================================
+  /*!
+   * \brief Construct the element box
+   */
+  //================================================================================
+
+  ElementBndBoxTree::ElementBox::ElementBox(const SMDS_MeshElement* elem)
+  {
+    _element  = elem;
+    _refCount = 1;
+    SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
+    while ( nIt->more() )
+      Add( TNodeXYZ( cast2Node( nIt->next() )));
+    Enlarge( NodeRadius );
+  }
+
+} // namespace
+
+//=======================================================================
+/*!
+ * \brief Implementation of search for the elements by point
+ */
+//=======================================================================
+
+struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
+{
+  SMESHDS_Mesh*           _mesh;
+  ElementBndBoxTree*      _ebbTree;
+  SMESH_NodeSearcherImpl* _nodeSearcher;
+  SMDSAbs_ElementType     _elementType;
+
+  SMESH_ElementSearcherImpl( SMESHDS_Mesh& mesh ): _mesh(&mesh),_ebbTree(0),_nodeSearcher(0) {}
+  ~SMESH_ElementSearcherImpl()
+  {
+    if ( _ebbTree )      delete _ebbTree;      _ebbTree      = 0;
+    if ( _nodeSearcher ) delete _nodeSearcher; _nodeSearcher = 0;
+  }
+
+  /*!
+   * \brief Return elements of given type where the given point is IN or ON.
+   *
+   * 'ALL' type means elements of any type excluding nodes and 0D elements
+   */
+  void FindElementsByPoint(const gp_Pnt&                      point,
+                           SMDSAbs_ElementType                type,
+                           vector< const SMDS_MeshElement* >& foundElements)
+  {
+    foundElements.clear();
+
+    const SMDS_MeshInfo& meshInfo = _mesh->GetMeshInfo();
+
+    // -----------------
+    // define tolerance
+    // -----------------
+    double tolerance = 0;
+    if ( _nodeSearcher && meshInfo.NbNodes() > 1 )
+    {
+      double boxSize = _nodeSearcher->getTree()->maxSize();
+      tolerance = 1e-8 * boxSize/* / meshInfo.NbNodes()*/;
+    }
+    else if ( _ebbTree && meshInfo.NbElements() > 0 )
+    {
+      double boxSize = _ebbTree->maxSize();
+      tolerance = 1e-8 * boxSize/* / meshInfo.NbElements()*/;
+    }
+    if ( tolerance == 0 )
+    {
+      // define tolerance by size of a most complex element
+      int complexType = SMDSAbs_Volume;
+      while ( complexType > SMDSAbs_All &&
+              meshInfo.NbElements( SMDSAbs_ElementType( complexType )) < 1 )
+        --complexType;
+      if ( complexType == SMDSAbs_All ) return; // empty mesh
+
+      double elemSize;
+      if ( complexType == int( SMDSAbs_Node ))
+      {
+        SMDS_NodeIteratorPtr nodeIt = _mesh->nodesIterator();
+        elemSize = 1;
+        if ( meshInfo.NbNodes() > 2 )
+          elemSize = TNodeXYZ( nodeIt->next() ).Distance( nodeIt->next() );
+      }
+      else
+      {
+        const SMDS_MeshElement* elem =
+          _mesh->elementsIterator( SMDSAbs_ElementType( complexType ))->next();
+        SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
+        TNodeXYZ n1( cast2Node( nodeIt->next() ));
+        while ( nodeIt->more() )
+        {
+          double dist = n1.Distance( cast2Node( nodeIt->next() ));
+          elemSize = max( dist, elemSize );
+        }
+      }
+      tolerance = 1e-6 * elemSize;
+    }
+
+    // =================================================================================
+    if ( type == SMDSAbs_Node || type == SMDSAbs_0DElement )
+    {
+      if ( !_nodeSearcher )
+        _nodeSearcher = new SMESH_NodeSearcherImpl( _mesh );
+
+      const SMDS_MeshNode* closeNode = _nodeSearcher->FindClosestTo( point );
+      if ( !closeNode ) return;
+
+      if ( point.Distance( TNodeXYZ( closeNode )) > tolerance )
+        return; // to far from any node
+
+      if ( type == SMDSAbs_Node )
+      {
+        foundElements.push_back( closeNode );
+      }
+      else
+      {
+        SMDS_ElemIteratorPtr elemIt = closeNode->GetInverseElementIterator( SMDSAbs_0DElement );
+        while ( elemIt->more() )
+          foundElements.push_back( elemIt->next() );
+      }
+    }
+    // =================================================================================
+    else // elements more complex than 0D
+    {
+      if ( !_ebbTree || _elementType != type )
+      {
+        if ( _ebbTree ) delete _ebbTree;
+        _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type );
+      }
+      TIDSortedElemSet suspectElems;
+      _ebbTree->getElementsNearPoint( point, suspectElems );
+      TIDSortedElemSet::iterator elem = suspectElems.begin();
+      for ( ; elem != suspectElems.end(); ++elem )
+        if ( !SMESH_MeshEditor::isOut( *elem, point, tolerance ))
+          foundElements.push_back( *elem );
+    }
+  }
+}; // struct SMESH_ElementSearcherImpl
+
+//=======================================================================
+/*!
+ * \brief Return SMESH_ElementSearcher
+ */
+//=======================================================================
+
+SMESH_ElementSearcher* SMESH_MeshEditor::GetElementSearcher()
+{
+  return new SMESH_ElementSearcherImpl( *GetMeshDS() );
+}
+
+//=======================================================================
+/*!
+ * \brief Return true if the point is IN or ON of the element
+ */
+//=======================================================================
+
+bool SMESH_MeshEditor::isOut( const SMDS_MeshElement* element, const gp_Pnt& point, double tol )
+{
+  if ( element->GetType() == SMDSAbs_Volume)
+  {
+    return SMDS_VolumeTool( element ).IsOut( point.X(), point.Y(), point.Z(), tol );
+  }
+
+  // get ordered nodes
+
+  vector< gp_XYZ > xyz;
+
+  SMDS_ElemIteratorPtr nodeIt = element->nodesIterator();
+  if ( element->IsQuadratic() )
+    if (const SMDS_QuadraticFaceOfNodes* f=dynamic_cast<const SMDS_QuadraticFaceOfNodes*>(element))
+      nodeIt = f->interlacedNodesElemIterator();
+    else if (const SMDS_QuadraticEdge*  e =dynamic_cast<const SMDS_QuadraticEdge*>(element))
+      nodeIt = e->interlacedNodesElemIterator();
+
+  while ( nodeIt->more() )
+    xyz.push_back( TNodeXYZ( cast2Node( nodeIt->next() )));
+
+  if ( element->GetType() == SMDSAbs_Face ) // --------------------------------------------------
+  {
+    // gravity center
+    gp_XYZ gc(0,0,0);
+    gc = accumulate( xyz.begin(), xyz.end(), gc );
+    gc /= element->NbNodes();
+
+    // compute face normal using gc
+    gp_Vec normal(0,0,0);
+    xyz.push_back( xyz.front() );
+    for ( int i = 0; i < element->NbNodes(); ++i )
+    {
+      gp_Vec edge( xyz[i], xyz[i+1]);
+      gp_Vec n2gc( xyz[i], gc );
+      normal += edge ^ n2gc;
+    }
+    double faceDoubleArea = normal.Magnitude();
+    if ( faceDoubleArea <= numeric_limits<double>::min() )
+      return true; // invalid face
+    normal /= faceDoubleArea;
+
+    // check if the point lays on face plane
+    gp_Vec n2p( xyz[0], point );
+    if ( fabs( n2p * normal ) > tol )
+      return true; // not on face plane
+
+    // check if point is out of face boundary
+    int i, out = false;
+    for ( i = 0; !out && i < element->NbNodes(); ++i )
+    {
+      gp_Vec edge( xyz[i], xyz[i+1]);
+      gp_Vec n2p ( xyz[i], point );
+      gp_Vec cross = edge ^ n2p;
+      out = ( cross * normal < -tol );
+    }
+    if ( out && element->IsPoly() )
+    {
+      // define point position by the closest edge
+      double minDist = numeric_limits<double>::max();
+      int iMinDist;
+      for ( i = 0; i < element->NbNodes(); ++i )
+      {
+        gp_Vec edge( xyz[i], xyz[i+1]);
+        gp_Vec n1p ( xyz[i], point);
+        double dist = ( edge ^ n1p ).Magnitude() / edge.Magnitude();
+        if ( dist < minDist )
+          iMinDist = i;
+      }
+      gp_Vec edge( xyz[iMinDist], xyz[iMinDist+1]);
+      gp_Vec n2p ( xyz[iMinDist], point );
+      gp_Vec cross = edge ^ n2p;
+      out = ( cross * normal < -tol );
+    }
+    return out;
+  }
+  if ( element->GetType() == SMDSAbs_Edge ) // --------------------------------------------------
+  {
+    for ( int i = 1; i < element->NbNodes(); ++i )
+    {
+      gp_Vec edge( xyz[i-1], xyz[i]);
+      gp_Vec n1p ( xyz[i-1], point);
+      double dist = ( edge ^ n1p ).Magnitude() / edge.Magnitude();
+      if ( dist > tol )
+        return true;
+      gp_Vec n2p( xyz[i], point );
+      if ( fabs( edge.Magnitude() - n1p.Magnitude() - n2p.Magnitude()) > tol )
+        return true;
+    }
+    return false;
+  }
+  // Node or 0D element -------------------------------------------------------------------------
+  {
+    gp_Vec n2p ( xyz[0], point );
+    return n2p.Magnitude() <= tol;
+  }
+  return true;
+}
+
 //=======================================================================
 //function : SimplifyFace
 //purpose  :
@@ -4818,7 +5777,7 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
       SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
       while ( invElemIt->more() ) {
         const SMDS_MeshElement* elem = invElemIt->next();
-          elems.insert(elem);
+        elems.insert(elem);
       }
     }
   }
@@ -4844,6 +5803,25 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
       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
         iRepl[ nbRepl++ ] = iCur;
       }
       curNodes[ iCur ] = n;
@@ -5333,24 +6311,24 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
 // ========================================================
 class SortableElement : public set <const SMDS_MeshElement*>
 {
- public:
+public:
 
   SortableElement( const SMDS_MeshElement* theElem )
-    {
-      myElem = theElem;
-      SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
-      while ( nodeIt->more() )
-        this->insert( nodeIt->next() );
-    }
+  {
+    myElem = theElem;
+    SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
+    while ( nodeIt->more() )
+      this->insert( nodeIt->next() );
+  }
 
   const SMDS_MeshElement* Get() const
-    { return myElem; }
+  { return myElem; }
 
   void Set(const SMDS_MeshElement* e) const
-    { myElem = e; }
+  { myElem = e; }
 
 
- private:
+private:
   mutable const SMDS_MeshElement* myElem;
 };
 
@@ -5360,7 +6338,7 @@ class SortableElement : public set <const SMDS_MeshElement*>
 //           Search among theElements or in the whole mesh if theElements is empty
 //=======================================================================
 void SMESH_MeshEditor::FindEqualElements(set<const SMDS_MeshElement*> & theElements,
-                                        TListOfListOfElementsID &      theGroupsOfElementsID)
+                                         TListOfListOfElementsID &      theGroupsOfElementsID)
 {
   myLastCreatedElems.Clear();
   myLastCreatedNodes.Clear();
@@ -5458,7 +6436,7 @@ void SMESH_MeshEditor::MergeElements(TListOfListOfElementsID & theGroupsOfElemen
 void SMESH_MeshEditor::MergeEqualElements()
 {
   set<const SMDS_MeshElement*> aMeshElements; /* empty input -
-                                                to merge equal elements in the whole mesh */
+                                                 to merge equal elements in the whole mesh */
   TListOfListOfElementsID aGroupsOfElementsID;
   FindEqualElements(aMeshElements, aGroupsOfElementsID);
   MergeElements(aGroupsOfElementsID);
@@ -5472,10 +6450,10 @@ void SMESH_MeshEditor::MergeEqualElements()
 //=======================================================================
 
 const SMDS_MeshElement*
-  SMESH_MeshEditor::FindFaceInSet(const SMDS_MeshNode*    n1,
-                                  const SMDS_MeshNode*    n2,
-                                  const TIDSortedElemSet& elemSet,
-                                  const TIDSortedElemSet& avoidSet)
+SMESH_MeshEditor::FindFaceInSet(const SMDS_MeshNode*    n1,
+                                const SMDS_MeshNode*    n2,
+                                const TIDSortedElemSet& elemSet,
+                                const TIDSortedElemSet& avoidSet)
 
 {
   SMDS_ElemIteratorPtr invElemIt = n1->GetInverseElementIterator(SMDSAbs_Face);
@@ -5589,7 +6567,7 @@ bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode*             theFirst
 
   //vector<const SMDS_MeshNode*> nodes;
   const SMDS_MeshNode *nIgnore = theFirstNode, *nStart = theSecondNode;
-  set < const SMDS_MeshElement* > foundElems;
+  TIDSortedElemSet foundElems;
   bool needTheLast = ( theLastNode != 0 );
 
   while ( nStart != theLastNode ) {
@@ -5608,7 +6586,7 @@ bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode*             theFirst
         int iNode = 0, nbNodes = e->NbNodes();
         //const SMDS_MeshNode* nodes[nbNodes+1];
         vector<const SMDS_MeshNode*> nodes(nbNodes+1);
-        
+
         if(e->IsQuadratic()) {
           const SMDS_QuadraticFaceOfNodes* F =
             static_cast<const SMDS_QuadraticFaceOfNodes*>(e);
@@ -5728,15 +6706,15 @@ bool SMESH_MeshEditor::CheckFreeBorderNodes(const SMDS_MeshNode* theNode1,
 //=======================================================================
 
 SMESH_MeshEditor::Sew_Error
-  SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
-                                   const SMDS_MeshNode* theBordSecondNode,
-                                   const SMDS_MeshNode* theBordLastNode,
-                                   const SMDS_MeshNode* theSideFirstNode,
-                                   const SMDS_MeshNode* theSideSecondNode,
-                                   const SMDS_MeshNode* theSideThirdNode,
-                                   const bool           theSideIsFreeBorder,
-                                   const bool           toCreatePolygons,
-                                   const bool           toCreatePolyedrs)
+SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
+                                 const SMDS_MeshNode* theBordSecondNode,
+                                 const SMDS_MeshNode* theBordLastNode,
+                                 const SMDS_MeshNode* theSideFirstNode,
+                                 const SMDS_MeshNode* theSideSecondNode,
+                                 const SMDS_MeshNode* theSideThirdNode,
+                                 const bool           theSideIsFreeBorder,
+                                 const bool           toCreatePolygons,
+                                 const bool           toCreatePolyedrs)
 {
   myLastCreatedElems.Clear();
   myLastCreatedNodes.Clear();
@@ -5990,7 +6968,7 @@ SMESH_MeshEditor::Sew_Error
 
   TListOfListOfNodes nodeGroupsToMerge;
   if ( nbNodes[0] == nbNodes[1] ||
-      ( theSideIsFreeBorder && !theSideThirdNode)) {
+       ( theSideIsFreeBorder && !theSideThirdNode)) {
 
     // all nodes are to be merged
 
@@ -6610,11 +7588,13 @@ void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode*        theBetweenNode
 //=======================================================================
 
 int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
-                                              SMESH_MesherHelper& theHelper,
-                                             const bool          theForce3d)
+                                             SMESH_MesherHelper& theHelper,
+                                             const bool          theForce3d)
 {
   int nbElem = 0;
   if( !theSm ) return nbElem;
+
+  const bool notFromGroups = false;
   SMDS_ElemIteratorPtr ElemItr = theSm->GetElements();
   while(ElemItr->more())
   {
@@ -6632,62 +7612,57 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
     }
     SMDSAbs_ElementType aType = elem->GetType();
 
-    theSm->RemoveElement(elem);
-    GetMeshDS()->SMDS_Mesh::RemoveFreeElement(elem);
+    GetMeshDS()->RemoveFreeElement(elem, theSm, notFromGroups);
 
     const SMDS_MeshElement* NewElem = 0;
 
     switch( aType )
     {
     case SMDSAbs_Edge :
-    {
-      NewElem = theHelper.AddEdge(aNds[0], aNds[1], id, theForce3d);
-      break;
-    }
+      {
+        NewElem = theHelper.AddEdge(aNds[0], aNds[1], id, theForce3d);
+        break;
+      }
     case SMDSAbs_Face :
-    {
-      switch(nbNodes)
       {
-      case 3:
-       NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
-       break;
-      case 4:
-       NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
-       break;
-      default:
-       continue;
+        switch(nbNodes)
+        {
+        case 3:
+          NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
+          break;
+        case 4:
+          NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
+          break;
+        default:
+          continue;
+        }
+        break;
       }
-      break;
-    }
     case SMDSAbs_Volume :
-    {
-      switch(nbNodes)
       {
-      case 4:
-       NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], id, true);
-       break;
-      case 6:
-       NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], aNds[4], aNds[5], id, true);
-       break;
-      case 8:
-       NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
-                                      aNds[4], aNds[5], aNds[6], aNds[7], id, true);
-       break;
-      default:
-       continue;
+        switch(nbNodes)
+        {
+        case 4:
+          NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
+          break;
+        case 6:
+          NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], aNds[4], aNds[5], id, theForce3d);
+          break;
+        case 8:
+          NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
+                                        aNds[4], aNds[5], aNds[6], aNds[7], id, theForce3d);
+          break;
+        default:
+          continue;
+        }
+        break;
       }
-      break;
-    }
     default :
       continue;
     }
+    ReplaceElemInGroups( elem, NewElem, GetMeshDS());
     if( NewElem )
-    {
-      AddToSameGroups( NewElem, elem, GetMeshDS());
       theSm->AddElement( NewElem );
-    }
-    if ( NewElem != elem )
-      RemoveElemFromGroups (elem, GetMeshDS());
   }
   return nbElem;
 }
@@ -6702,6 +7677,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
 
   SMESH_MesherHelper aHelper(*myMesh);
   aHelper.SetIsQuadratic( true );
+  const bool notFromGroups = false;
 
   int nbCheckedElems = 0;
   if ( myMesh->HasShapeToMesh() )
@@ -6713,31 +7689,30 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
         SMESH_subMesh* sm = smIt->next();
         if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() ) {
           aHelper.SetSubShape( sm->GetSubShape() );
+          if ( !theForce3d) aHelper.SetCheckNodePosition(true);
           nbCheckedElems += convertElemToQuadratic(smDS, aHelper, theForce3d);
         }
       }
     }
   }
   int totalNbElems = meshDS->NbEdges() + meshDS->NbFaces() + meshDS->NbVolumes();
-  if ( nbCheckedElems < totalNbElems ) // not all elements in submeshes
+  if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes
   {
+    SMESHDS_SubMesh *smDS = 0;
     SMDS_EdgeIteratorPtr aEdgeItr = meshDS->edgesIterator();
     while(aEdgeItr->more())
     {
       const SMDS_MeshEdge* edge = aEdgeItr->next();
       if(edge && !edge->IsQuadratic())
       {
-       int id = edge->GetID();
-       const SMDS_MeshNode* n1 = edge->GetNode(0);
-       const SMDS_MeshNode* n2 = edge->GetNode(1);
+        int id = edge->GetID();
+        const SMDS_MeshNode* n1 = edge->GetNode(0);
+        const SMDS_MeshNode* n2 = edge->GetNode(1);
 
-       meshDS->SMDS_Mesh::RemoveFreeElement(edge);
+        meshDS->RemoveFreeElement(edge, smDS, notFromGroups);
 
         const SMDS_MeshEdge* NewEdge = aHelper.AddEdge(n1, n2, id, theForce3d);
-        if ( NewEdge )
-          AddToSameGroups(NewEdge, edge, meshDS);
-        if ( NewEdge != edge )
-          RemoveElemFromGroups (edge, meshDS);
+        ReplaceElemInGroups( edge, NewEdge, GetMeshDS());
       }
     }
     SMDS_FaceIteratorPtr aFaceItr = meshDS->facesIterator();
@@ -6752,27 +7727,24 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
 
       for(int i = 0; i < nbNodes; i++)
       {
-       aNds[i] = face->GetNode(i);
+        aNds[i] = face->GetNode(i);
       }
 
-      meshDS->SMDS_Mesh::RemoveFreeElement(face);
+      meshDS->RemoveFreeElement(face, smDS, notFromGroups);
 
       SMDS_MeshFace * NewFace = 0;
       switch(nbNodes)
       {
       case 3:
-       NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
-       break;
+        NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
+        break;
       case 4:
-       NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
-       break;
+        NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
+        break;
       default:
-       continue;
+        continue;
       }
-      if ( NewFace )
-        AddToSameGroups(NewFace, face, meshDS);
-      if ( NewFace != face )
-        RemoveElemFromGroups (face, meshDS);
+      ReplaceElemInGroups( face, NewFace, GetMeshDS());
     }
     SMDS_VolumeIteratorPtr aVolumeItr = meshDS->volumesIterator();
     while(aVolumeItr->more())
@@ -6786,35 +7758,36 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
 
       for(int i = 0; i < nbNodes; i++)
       {
-       aNds[i] = volume->GetNode(i);
+        aNds[i] = volume->GetNode(i);
       }
 
-      meshDS->SMDS_Mesh::RemoveFreeElement(volume);
+      meshDS->RemoveFreeElement(volume, smDS, notFromGroups);
 
       SMDS_MeshVolume * NewVolume = 0;
       switch(nbNodes)
       {
       case 4:
-       NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
-                                      aNds[3], id, true );
-       break;
+        NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
+                                      aNds[3], id, theForce3d );
+        break;
       case 6:
-       NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
-                                      aNds[3], aNds[4], aNds[5], id, true);
-       break;
+        NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
+                                      aNds[3], aNds[4], aNds[5], id, theForce3d);
+        break;
       case 8:
-       NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
-                                      aNds[4], aNds[5], aNds[6], aNds[7], id, true);
-       break;
+        NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
+                                      aNds[4], aNds[5], aNds[6], aNds[7], id, theForce3d);
+        break;
       default:
-       continue;
+        continue;
       }
-      if ( NewVolume )
-        AddToSameGroups(NewVolume, volume, meshDS);
-      if ( NewVolume != volume )
-        RemoveElemFromGroups (volume, meshDS);
+      ReplaceElemInGroups(volume, NewVolume, meshDS);
     }
   }
+  if ( !theForce3d ) {
+    aHelper.SetSubShape(0); // apply to the whole mesh
+    aHelper.FixQuadraticElements();
+  }
 }
 
 //=======================================================================
@@ -6830,6 +7803,8 @@ int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh *    theSm,
 {
   int nbElem = 0;
   SMESHDS_Mesh* meshDS = GetMeshDS();
+  const bool notFromGroups = false;
+
   while( theItr->more() )
   {
     const SMDS_MeshElement* elem = theItr->next();
@@ -6844,40 +7819,35 @@ int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh *    theSm,
 
       for(int i = 0; i < nbNodes; i++)
       {
-       const SMDS_MeshNode* n = elem->GetNode(i);
+        const SMDS_MeshNode* n = elem->GetNode(i);
 
-       if( elem->IsMediumNode( n ) )
+        if( elem->IsMediumNode( n ) )
           mediumNodes.push_back( n );
-       else
-         aNds.push_back( n );
+        else
+          aNds.push_back( n );
       }
       if( aNds.empty() ) continue;
       SMDSAbs_ElementType aType = elem->GetType();
 
       //remove old quadratic element
-      meshDS->SMDS_Mesh::RemoveFreeElement( elem );
-      if ( theSm )
-        theSm->RemoveElement( elem );
+      meshDS->RemoveFreeElement( elem, theSm, notFromGroups );
 
       SMDS_MeshElement * NewElem = AddElement( aNds, aType, false, id );
-      if ( NewElem )
-        AddToSameGroups(NewElem, elem, meshDS);
-      if ( NewElem != elem )
-        RemoveElemFromGroups (elem, meshDS);
+      ReplaceElemInGroups(elem, NewElem, meshDS);
       if( theSm && NewElem )
-       theSm->AddElement( NewElem );
+        theSm->AddElement( NewElem );
 
       // remove medium nodes
       vector<const SMDS_MeshNode*>::iterator nIt = mediumNodes.begin();
       for ( ; nIt != mediumNodes.end(); ++nIt ) {
         const SMDS_MeshNode* n = *nIt;
-        if ( n->NbInverseNodes() == 0 ) {
+        if ( n->NbInverseElements() == 0 ) {
           if ( n->GetPosition()->GetShapeId() != theShapeID )
             meshDS->RemoveFreeNode( n, meshDS->MeshElements
                                     ( n->GetPosition()->GetShapeId() ));
           else
             meshDS->RemoveFreeNode( n, theSm );
-       }
+        }
       }
     }
   }
@@ -6903,10 +7873,10 @@ bool  SMESH_MeshEditor::ConvertFromQuadratic()
       }
     }
   }
-  
+
   int totalNbElems =
     GetMeshDS()->NbEdges() + GetMeshDS()->NbFaces() + GetMeshDS()->NbVolumes();
-  if ( nbCheckedElems < totalNbElems ) // not all elements in submeshes
+  if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes
   {
     SMESHDS_SubMesh *aSM = 0;
     removeQuadElem( aSM, GetMeshDS()->elementsIterator(), 0 );
@@ -6921,12 +7891,12 @@ bool  SMESH_MeshEditor::ConvertFromQuadratic()
 //=======================================================================
 
 SMESH_MeshEditor::Sew_Error
-  SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
-                                     TIDSortedElemSet&    theSide2,
-                                     const SMDS_MeshNode* theFirstNode1,
-                                     const SMDS_MeshNode* theFirstNode2,
-                                     const SMDS_MeshNode* theSecondNode1,
-                                     const SMDS_MeshNode* theSecondNode2)
+SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
+                                   TIDSortedElemSet&    theSide2,
+                                   const SMDS_MeshNode* theFirstNode1,
+                                   const SMDS_MeshNode* theFirstNode2,
+                                   const SMDS_MeshNode* theSecondNode1,
+                                   const SMDS_MeshNode* theSecondNode2)
 {
   myLastCreatedElems.Clear();
   myLastCreatedNodes.Clear();
@@ -7165,40 +8135,40 @@ SMESH_MeshEditor::Sew_Error
           const SMDS_MeshElement* aFreeFace = freeFaceList.front();
           faceSet->insert( aFreeFace );
           // complete a node set with nodes of a found free face
-//           for ( iNode = 0; iNode < ; iNode++ )
-//             nodeSet->insert( fNodes[ iNode ] );
+          //           for ( iNode = 0; iNode < ; iNode++ )
+          //             nodeSet->insert( fNodes[ iNode ] );
         }
 
       } // loop on volumes of a side
 
-//       // complete a set of faces if new nodes in a nodeSet appeared
-//       // ----------------------------------------------------------
-//       if ( nodeSetSize != nodeSet->size() ) {
-//         for ( ; nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
-//           SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
-//           while ( fIt->more() ) { // loop on faces sharing a node
-//             const SMDS_MeshElement* f = fIt->next();
-//             if ( faceSet->find( f ) == faceSet->end() ) {
-//               // check if all nodes are in nodeSet and
-//               // complete setOfFaceNodeSet if they are
-//               set <const SMDS_MeshNode*> faceNodeSet;
-//               SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
-//               bool allInSet = true;
-//               while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
-//                 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
-//                 if ( nodeSet->find( n ) == nodeSet->end() )
-//                   allInSet = false;
-//                 else
-//                   faceNodeSet.insert( n );
-//               }
-//               if ( allInSet ) {
-//                 faceSet->insert( f );
-//                 setOfFaceNodeSet.insert( faceNodeSet );
-//               }
-//             }
-//           }
-//         }
-//       }
+      //       // complete a set of faces if new nodes in a nodeSet appeared
+      //       // ----------------------------------------------------------
+      //       if ( nodeSetSize != nodeSet->size() ) {
+      //         for ( ; nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
+      //           SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
+      //           while ( fIt->more() ) { // loop on faces sharing a node
+      //             const SMDS_MeshElement* f = fIt->next();
+      //             if ( faceSet->find( f ) == faceSet->end() ) {
+      //               // check if all nodes are in nodeSet and
+      //               // complete setOfFaceNodeSet if they are
+      //               set <const SMDS_MeshNode*> faceNodeSet;
+      //               SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
+      //               bool allInSet = true;
+      //               while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
+      //                 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+      //                 if ( nodeSet->find( n ) == nodeSet->end() )
+      //                   allInSet = false;
+      //                 else
+      //                   faceNodeSet.insert( n );
+      //               }
+      //               if ( allInSet ) {
+      //                 faceSet->insert( f );
+      //                 setOfFaceNodeSet.insert( faceNodeSet );
+      //               }
+      //             }
+      //           }
+      //         }
+      //       }
     } // Create temporary faces, if there are volumes given
   } // loop on sides
 
@@ -7304,10 +8274,10 @@ SMESH_MeshEditor::Sew_Error
                   nbl++;
                   if(iSide==0)
                     notLinkNodes1[nbl] = n;
-                    //notLinkNodes1.push_back(n);
+                  //notLinkNodes1.push_back(n);
                   else
                     notLinkNodes2[nbl] = n;
-                    //notLinkNodes2.push_back(n);
+                  //notLinkNodes2.push_back(n);
                 }
                 //faceNodes[ iSide ][ iNode++ ] = n;
                 if(iSide==0) {
@@ -7388,7 +8358,7 @@ SMESH_MeshEditor::Sew_Error
         //nReplaceMap.insert( TNodeNodeMap::value_type
         //                   ( notLinkNodes[0][0], notLinkNodes[1][0] ));
         nReplaceMap.insert( TNodeNodeMap::value_type
-                           ( notLinkNodes1[0], notLinkNodes2[0] ));
+                            ( notLinkNodes1[0], notLinkNodes2[0] ));
       }
       else {
         for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
@@ -7409,7 +8379,7 @@ SMESH_MeshEditor::Sew_Error
           //                   ( notLinkNodes[0][1], notLinkNodes[1][1] ));
           for(int nn=0; nn<nbNodes-2; nn++) {
             nReplaceMap.insert( TNodeNodeMap::value_type
-                             ( notLinkNodes1[nn], notLinkNodes2[nn] ));
+                                ( notLinkNodes1[nn], notLinkNodes2[nn] ));
           }
         }
         else {
@@ -7419,7 +8389,7 @@ SMESH_MeshEditor::Sew_Error
           //                   ( notLinkNodes[0][1], notLinkNodes[1][0] ));
           for(int nn=0; nn<nbNodes-2; nn++) {
             nReplaceMap.insert( TNodeNodeMap::value_type
-                             ( notLinkNodes1[nn], notLinkNodes2[nbNodes-3-nn] ));
+                                ( notLinkNodes1[nn], notLinkNodes2[nbNodes-3-nn] ));
           }
         }
       }
@@ -7449,10 +8419,10 @@ SMESH_MeshEditor::Sew_Error
   } // loop on link lists
 
   if ( aResult == SEW_OK &&
-      ( linkIt[0] != linkList[0].end() ||
-       !faceSetPtr[0]->empty() || !faceSetPtr[1]->empty() )) {
+       ( linkIt[0] != linkList[0].end() ||
+         !faceSetPtr[0]->empty() || !faceSetPtr[1]->empty() )) {
     MESSAGE( (linkIt[0] != linkList[0].end()) <<" "<< (faceSetPtr[0]->empty()) <<
-            " " << (faceSetPtr[1]->empty()));
+             " " << (faceSetPtr[1]->empty()));
     aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
   }
 
@@ -7507,17 +8477,17 @@ SMESH_MeshEditor::Sew_Error
 }
 
 //================================================================================
-  /*!
  * \brief Find corresponding nodes in two sets of faces
   * \param theSide1 - first face set
   * \param theSide2 - second first face
   * \param theFirstNode1 - a boundary node of set 1
   * \param theFirstNode2 - a node of set 2 corresponding to theFirstNode1
   * \param theSecondNode1 - a boundary node of set 1 linked with theFirstNode1
   * \param theSecondNode2 - a node of set 2 corresponding to theSecondNode1
   * \param nReplaceMap - output map of corresponding nodes
   * \retval bool  - is a success or not
  */
+/*!
+ * \brief Find corresponding nodes in two sets of faces
+ * \param theSide1 - first face set
+ * \param theSide2 - second first face
+ * \param theFirstNode1 - a boundary node of set 1
+ * \param theFirstNode2 - a node of set 2 corresponding to theFirstNode1
+ * \param theSecondNode1 - a boundary node of set 1 linked with theFirstNode1
+ * \param theSecondNode2 - a node of set 2 corresponding to theSecondNode1
+ * \param nReplaceMap - output map of corresponding nodes
+ * \retval bool  - is a success or not
+ */
 //================================================================================
 
 #ifdef _DEBUG_
@@ -7541,8 +8511,8 @@ SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
   if ( theSecondNode1 != theSecondNode2 )
     nReplaceMap.insert( make_pair( theSecondNode1, theSecondNode2 ));
 
-  set< TLink > linkSet; // set of nodes where order of nodes is ignored
-  linkSet.insert( TLink( theFirstNode1, theSecondNode1 ));
+  set< SMESH_TLink > linkSet; // set of nodes where order of nodes is ignored
+  linkSet.insert( SMESH_TLink( theFirstNode1, theSecondNode1 ));
 
   list< NLink > linkList[2];
   linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
@@ -7633,10 +8603,9 @@ SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
         return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
       }
 #ifdef DEBUG_MATCHING_NODES
-      cout << " Link 1: " << link[0].first->GetID() <<" "<< link[0].second->GetID()
-           << " F 1: " << face[0];
-      cout << "| Link 2: " << link[1].first->GetID() <<" "<< link[1].second->GetID()
-           << " F 2: " << face[1] << " | Bind: "<<endl ;
+      MESSAGE ( " Link 1: " << link[0].first->GetID() <<" "<< link[0].second->GetID()
+                << " F 1: " << face[0] << "| Link 2: " << link[1].first->GetID() <<" "
+                << link[1].second->GetID() << " F 2: " << face[1] << " | Bind: " ) ;
 #endif
       int nbN = nbNodes[0];
       {
@@ -7644,7 +8613,7 @@ SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
         list<const SMDS_MeshNode*>::iterator n2 = notLinkNodes[1].begin();
         for ( int i = 0 ; i < nbN - 2; ++i ) {
 #ifdef DEBUG_MATCHING_NODES
-          cout << (*n1)->GetID() << " to " << (*n2)->GetID() << endl;
+          MESSAGE ( (*n1)->GetID() << " to " << (*n2)->GetID() );
 #endif
           nReplaceMap.insert( make_pair( *(n1++), *(n2++) ));
         }
@@ -7658,16 +8627,16 @@ SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
       for ( int i = 0; i < nbN; i++ )
       {
         const SMDS_MeshNode* n2 = f0->GetNode( i );
-        pair< set< TLink >::iterator, bool > iter_isnew =
-          linkSet.insert( TLink( n1, n2 ));
+        pair< set< SMESH_TLink >::iterator, bool > iter_isnew =
+          linkSet.insert( SMESH_TLink( n1, n2 ));
         if ( !iter_isnew.second ) { // already in a set: no need to process
           linkSet.erase( iter_isnew.first );
         }
         else // new in set == encountered for the first time: add
         {
 #ifdef DEBUG_MATCHING_NODES
-          cout << "Add link 1: " << n1->GetID() << " " << n2->GetID() << " ";
-          cout << " | link 2: " << nReplaceMap[n1]->GetID() << " " << nReplaceMap[n2]->GetID() << " " << endl;
+          MESSAGE ( "Add link 1: " << n1->GetID() << " " << n2->GetID() << " "
+                    << " | link 2: " << nReplaceMap[n1]->GetID() << " " << nReplaceMap[n2]->GetID() << " " );
 #endif
           linkList[0].push_back ( NLink( n1, n2 ));
           linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
@@ -7679,3 +8648,222 @@ SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
 
   return SEW_OK;
 }
+
+/*!
+  \brief Creates a hole in a mesh by doubling the nodes of some particular elements
+  \param theElems - the list of elements (edges or faces) to be replicated
+  The nodes for duplication could be found from these elements
+  \param theNodesNot - list of nodes to NOT replicate
+  \param theAffectedElems - the list of elements (cells and edges) to which the 
+  replicated nodes should be associated to.
+  \return TRUE if operation has been completed successfully, FALSE otherwise
+*/
+bool SMESH_MeshEditor::DoubleNodes( const TIDSortedElemSet& theElems,
+                                    const TIDSortedElemSet& theNodesNot,
+                                    const TIDSortedElemSet& theAffectedElems )
+{
+  myLastCreatedElems.Clear();
+  myLastCreatedNodes.Clear();
+
+  if ( theElems.size() == 0 )
+    return false;
+
+  SMESHDS_Mesh* aMeshDS = GetMeshDS();
+  if ( !aMeshDS )
+    return false;
+
+  bool res = false;
+  std::map< const SMDS_MeshNode*, const SMDS_MeshNode* > anOldNodeToNewNode;
+  // duplicate elements and nodes
+  res = doubleNodes( aMeshDS, theElems, theNodesNot, anOldNodeToNewNode, true );
+  // replce nodes by duplications
+  res = doubleNodes( aMeshDS, theAffectedElems, theNodesNot, anOldNodeToNewNode, false );
+  return res;
+}
+
+/*!
+  \brief Creates a hole in a mesh by doubling the nodes of some particular elements
+  \param theMeshDS - mesh instance
+  \param theElems - the elements replicated or modified (nodes should be changed)
+  \param theNodesNot - nodes to NOT replicate
+  \param theNodeNodeMap - relation of old node to new created node
+  \param theIsDoubleElem - flag os to replicate element or modify
+  \return TRUE if operation has been completed successfully, FALSE otherwise
+*/
+bool SMESH_MeshEditor::doubleNodes( SMESHDS_Mesh*     theMeshDS,
+                                    const TIDSortedElemSet& theElems,
+                                    const TIDSortedElemSet& theNodesNot,
+                                    std::map< const SMDS_MeshNode*,
+                                    const SMDS_MeshNode* >& theNodeNodeMap,
+                                    const bool theIsDoubleElem )
+{
+  // iterate on through element and duplicate them (by nodes duplication)
+  bool res = false;
+  TIDSortedElemSet::const_iterator elemItr = theElems.begin();
+  for ( ;  elemItr != theElems.end(); ++elemItr )
+  {
+    const SMDS_MeshElement* anElem = *elemItr;
+    if (!anElem)
+      continue;
+
+    bool isDuplicate = false;
+    // duplicate nodes to duplicate element
+    std::vector<const SMDS_MeshNode*> newNodes( anElem->NbNodes() );
+    SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
+    int ind = 0;
+    while ( anIter->more() ) 
+    { 
+
+      SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next();
+      SMDS_MeshNode* aNewNode = aCurrNode;
+      if ( theNodeNodeMap.find( aCurrNode ) != theNodeNodeMap.end() )
+        aNewNode = (SMDS_MeshNode*)theNodeNodeMap[ aCurrNode ];
+      else if ( theIsDoubleElem && theNodesNot.find( aCurrNode ) == theNodesNot.end() )
+      {
+        // duplicate node
+        aNewNode = theMeshDS->AddNode( aCurrNode->X(), aCurrNode->Y(), aCurrNode->Z() );
+        theNodeNodeMap[ aCurrNode ] = aNewNode;
+        myLastCreatedNodes.Append( aNewNode );
+      }
+      isDuplicate |= (aCurrNode == aNewNode);
+      newNodes[ ind++ ] = aNewNode;
+    }
+    if ( !isDuplicate )
+      continue;
+
+    if ( theIsDoubleElem )
+      myLastCreatedElems.Append( AddElement(newNodes, anElem->GetType(), anElem->IsPoly()) );
+    else
+      theMeshDS->ChangeElementNodes( anElem, &newNodes[ 0 ], anElem->NbNodes() );
+
+    res = true;
+  }
+  return res;
+}
+
+/*!
+  \brief Check if element located inside shape
+  \return TRUE if IN or ON shape, FALSE otherwise
+*/
+
+static bool isInside(const SMDS_MeshElement* theElem,
+                     BRepClass3d_SolidClassifier& theBsc3d,
+                     const double theTol)
+{
+  gp_XYZ centerXYZ (0, 0, 0);
+  SMDS_ElemIteratorPtr aNodeItr = theElem->nodesIterator();
+  while (aNodeItr->more())
+  {
+    SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
+    centerXYZ += gp_XYZ(aNode->X(), aNode->Y(), aNode->Z());
+  }
+  gp_Pnt aPnt(centerXYZ);
+  theBsc3d.Perform(aPnt, theTol);
+  TopAbs_State aState = theBsc3d.State();
+  return (aState == TopAbs_IN || aState == TopAbs_ON );
+}
+
+/*!
+  \brief Creates a hole in a mesh by doubling the nodes of some particular elements
+  \param theElems - group of of elements (edges or faces) to be replicated
+  \param theNodesNot - group of nodes not to replicated
+  \param theShape - shape to detect affected elements (element which geometric center
+  located on or inside shape).
+  The replicated nodes should be associated to affected elements.
+  \return TRUE if operation has been completed successfully, FALSE otherwise
+*/
+
+bool SMESH_MeshEditor::DoubleNodesInRegion( const TIDSortedElemSet& theElems,
+                                            const TIDSortedElemSet& theNodesNot,
+                                            const TopoDS_Shape&     theShape )
+{
+  if ( theShape.IsNull() )
+    return false;
+
+  const double aTol = Precision::Confusion();
+  BRepClass3d_SolidClassifier bsc3d(theShape);
+  bsc3d.PerformInfinitePoint(aTol);
+
+  // iterates on indicated elements and get elements by back references from their nodes
+  TIDSortedElemSet anAffected;
+  TIDSortedElemSet::const_iterator elemItr = theElems.begin();
+  for ( ;  elemItr != theElems.end(); ++elemItr )
+  {
+    SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr;
+    if (!anElem)
+      continue;
+
+    SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator();
+    while ( nodeItr->more() )
+    {
+      const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>(nodeItr->next());
+      if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() )
+        continue;
+      SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator();
+      while ( backElemItr->more() )
+      {
+        SMDS_MeshElement* curElem = (SMDS_MeshElement*)backElemItr->next();
+        if ( curElem && theElems.find(curElem) == theElems.end() &&
+             isInside( curElem, bsc3d, aTol ) )
+          anAffected.insert( curElem );
+      }
+    }
+  }
+  return DoubleNodes( theElems, theNodesNot, anAffected );
+}
+
+/*!
+ * \brief Generated skin mesh (containing 2D cells) from 3D mesh
+ * The created 2D mesh elements based on nodes of free faces of boundary volumes
+ * \return TRUE if operation has been completed successfully, FALSE otherwise
+ */
+
+bool SMESH_MeshEditor::Make2DMeshFrom3D()
+{
+  // iterates on volume elements and detect all free faces on them
+  SMESHDS_Mesh* aMesh = GetMeshDS();
+  if (!aMesh)
+    return false;
+  bool res = false;
+  SMDS_VolumeIteratorPtr vIt = aMesh->volumesIterator();
+  while(vIt->more())
+  {
+    const SMDS_MeshVolume* volume = vIt->next();
+    SMDS_VolumeTool vTool( volume );
+    const bool isPoly = volume->IsPoly();
+    const bool isQuad = volume->IsQuadratic();
+    for ( int iface = 0, n = vTool.NbFaces(); iface < n; iface++ )
+    {
+      if (!vTool.IsFreeFace(iface))
+        continue;
+      vector<const SMDS_MeshNode *> nodes;
+      int nbFaceNodes = vTool.NbFaceNodes(iface);
+      const SMDS_MeshNode** faceNodes = vTool.GetFaceNodes(iface);
+      if (vTool.IsFaceExternal(iface)) 
+      {
+        int inode = 0;
+        for ( ; inode < nbFaceNodes; inode += isQuad ? 2 : 1)
+          nodes.push_back(faceNodes[inode]);
+        if (isQuad)
+          for ( inode = 1; inode < nbFaceNodes; inode += 2)
+            nodes.push_back(faceNodes[inode]);
+      }
+      else
+      {
+        int inode = nbFaceNodes-1;
+        for ( ; inode >=0; inode -= isQuad ? 2 : 1)
+          nodes.push_back(faceNodes[inode]);
+        if (isQuad)
+          for ( inode = nbFaceNodes-2; inode >=0; inode -= 2)
+            nodes.push_back(faceNodes[inode]);
+      }
+
+      // add new face based on volume nodes
+      if (aMesh->FindFace( nodes ) )
+        continue; // face already exsist
+      myLastCreatedElems.Append( AddElement(nodes, SMDSAbs_Face, isPoly && iface == 1) );
+      res = true;
+    }
+  }
+  return res;
+}