Salome HOME
Copyright update 2022
[modules/smesh.git] / src / SMESH / SMESH_MeshEditor.cxx
index 81db6b65f296d5cd3ed12332810c2732bd98efcc..be874359b73c68b2ab7891c3acb3c09104a2984c 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2022  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
 
 #include "SMESH_TryCatch.hxx" // include after OCCT headers!
 
 
 #include "SMESH_TryCatch.hxx" // include after OCCT headers!
 
+#include <smIdType.hxx>
+
 #define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
 
 using namespace std;
 #define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
 
 using namespace std;
@@ -161,8 +163,7 @@ SMESH_MeshEditor::ElemFeatures::Init( const SMDS_MeshElement* elem, bool basicOn
         myIsQuad = elem->IsQuadratic();
         if ( myType == SMDSAbs_Volume && !basicOnly )
         {
         myIsQuad = elem->IsQuadratic();
         if ( myType == SMDSAbs_Volume && !basicOnly )
         {
-          vector<int> quant = static_cast<const SMDS_MeshVolume* >( elem )->GetQuantities();
-          myPolyhedQuantities.swap( quant );
+          myPolyhedQuantities = static_cast<const SMDS_MeshVolume* >( elem )->GetQuantities();
         }
       }
     }
         }
       }
     }
@@ -187,7 +188,7 @@ SMESH_MeshEditor::AddElement(const vector<const SMDS_MeshNode*> & node,
   SMDS_MeshElement* e = 0;
   int nbnode = node.size();
   SMESHDS_Mesh* mesh = GetMeshDS();
   SMDS_MeshElement* e = 0;
   int nbnode = node.size();
   SMESHDS_Mesh* mesh = GetMeshDS();
-  const int ID = features.myID;
+  const smIdType ID = features.myID;
 
   switch ( features.myType ) {
   case SMDSAbs_Face:
 
   switch ( features.myType ) {
   case SMDSAbs_Face:
@@ -297,6 +298,18 @@ SMESH_MeshEditor::AddElement(const vector<const SMDS_MeshNode*> & node,
                                                  node[8], node[9], node[10],node[11],
                                                  node[12],node[13],node[14] );
       }
                                                  node[8], node[9], node[10],node[11],
                                                  node[12],node[13],node[14] );
       }
+      else if (nbnode == 18) {
+        if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
+                                                 node[4], node[5], node[6], node[7],
+                                                 node[8], node[9], node[10],node[11],
+                                                 node[12],node[13],node[14],
+                                                 node[15],node[16],node[17],ID );
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
+                                                 node[4], node[5], node[6], node[7],
+                                                 node[8], node[9], node[10],node[11],
+                                                 node[12],node[13],node[14],
+                                                 node[15],node[16],node[17] );
+      }
       else if (nbnode == 20) {
         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
                                                  node[4], node[5], node[6], node[7],
       else if (nbnode == 20) {
         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
                                                  node[4], node[5], node[6], node[7],
@@ -378,12 +391,12 @@ SMESH_MeshEditor::AddElement(const vector<const SMDS_MeshNode*> & node,
  */
 //=======================================================================
 
  */
 //=======================================================================
 
-SMDS_MeshElement* SMESH_MeshEditor::AddElement(const vector<int> & nodeIDs,
-                                               const ElemFeatures& features)
+SMDS_MeshElement* SMESH_MeshEditor::AddElement(const vector<smIdType> & nodeIDs,
+                                               const ElemFeatures&      features)
 {
   vector<const SMDS_MeshNode*> nodes;
   nodes.reserve( nodeIDs.size() );
 {
   vector<const SMDS_MeshNode*> nodes;
   nodes.reserve( nodeIDs.size() );
-  vector<int>::const_iterator id = nodeIDs.begin();
+  vector<smIdType>::const_iterator id = nodeIDs.begin();
   while ( id != nodeIDs.end() ) {
     if ( const SMDS_MeshNode* node = GetMeshDS()->FindNode( *id++ ))
       nodes.push_back( node );
   while ( id != nodeIDs.end() ) {
     if ( const SMDS_MeshNode* node = GetMeshDS()->FindNode( *id++ ))
       nodes.push_back( node );
@@ -399,16 +412,16 @@ SMDS_MeshElement* SMESH_MeshEditor::AddElement(const vector<int> & nodeIDs,
 //           Modify a compute state of sub-meshes which become empty
 //=======================================================================
 
 //           Modify a compute state of sub-meshes which become empty
 //=======================================================================
 
-int SMESH_MeshEditor::Remove (const list< int >& theIDs,
-                              const bool         isNodes )
+smIdType SMESH_MeshEditor::Remove (const std::list< smIdType >& theIDs,
+                                   const bool                   isNodes )
 {
   ClearLastCreated();
 
   SMESHDS_Mesh* aMesh = GetMeshDS();
   set< SMESH_subMesh *> smmap;
 
 {
   ClearLastCreated();
 
   SMESHDS_Mesh* aMesh = GetMeshDS();
   set< SMESH_subMesh *> smmap;
 
-  int removed = 0;
-  list<int>::const_iterator it = theIDs.begin();
+  smIdType removed = 0;
+  list<smIdType>::const_iterator it = theIDs.begin();
   for ( ; it != theIDs.end(); it++ ) {
     const SMDS_MeshElement * elem;
     if ( isNodes )
   for ( ; it != theIDs.end(); it++ ) {
     const SMDS_MeshElement * elem;
     if ( isNodes )
@@ -461,10 +474,87 @@ int SMESH_MeshEditor::Remove (const list< int >& theIDs,
   return removed;
 }
 
   return removed;
 }
 
+//================================================================================
+/*!
+ * \brief Remove a node and fill a hole appeared, by changing surrounding faces
+ */
+//================================================================================
+
+void SMESH_MeshEditor::RemoveNodeWithReconnection( const SMDS_MeshNode* node )
+{
+  if ( ! node )
+    return;
+
+  if ( node->NbInverseElements( SMDSAbs_Volume ) > 0 )
+    throw SALOME_Exception( "RemoveNodeWithReconnection() applies to 2D mesh only" );
+
+  // check that only triangles surround the node
+  for ( SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator( SMDSAbs_Face ); fIt->more(); )
+  {
+    const SMDS_MeshElement* face = fIt->next();
+    if ( face->GetGeomType() != SMDSGeom_TRIANGLE )
+      throw SALOME_Exception( "RemoveNodeWithReconnection() applies to triangle mesh only" );
+    if ( face->IsQuadratic() )
+      throw SALOME_Exception( "RemoveNodeWithReconnection() applies to linear mesh only" );
+  }
+
+  std::vector< const SMDS_MeshNode*> neighbours(2);
+  SMESH_MeshAlgos::IsOn2DBoundary( node, & neighbours );
+
+  bool toRemove = ( neighbours.size() > 2 ); // non-manifold ==> just remove
+
+  // if ( neighbours.size() == 2 ) // on boundary
+  // {
+  //   // check if theNode and neighbours are on a line
+  //   gp_Pnt pN = SMESH_NodeXYZ( node );
+  //   gp_Pnt p0 = SMESH_NodeXYZ( neighbours[0] );
+  //   gp_Pnt p1 = SMESH_NodeXYZ( neighbours[1] );
+  //   double dist01 = p0.Distance( p1 );
+  //   double    tol = 0.01 * dist01;
+  //   double  distN = ( gp_Vec( p0, p1 ) ^ gp_Vec( p0, pN )).Magnitude() / dist01;
+  //   bool   onLine = distN < tol;
+  //   toRemove = !onLine;
+  // }
+
+  if ( neighbours.empty() ) // not on boundary
+  {
+    TIDSortedElemSet linkedNodes;
+    GetLinkedNodes( node, linkedNodes, SMDSAbs_Face );
+    for ( const SMDS_MeshElement* e : linkedNodes ) neighbours.push_back( cast2Node( e ));
+    if ( neighbours.empty() )
+      toRemove = true;
+  }
+
+  if ( toRemove )
+  {
+    this->Remove( std::list< smIdType >( 1, node->GetID() ), /*isNode=*/true );
+    return;
+  }
+
+  // choose a node to replace by
+  const SMDS_MeshNode* nToReplace = nullptr;
+  SMESH_NodeXYZ           nodeXYZ = node;
+  double                  minDist = Precision::Infinite();
+  for ( const SMDS_MeshNode* n : neighbours )
+  {
+    double dist = nodeXYZ.SquareDistance( n );
+    if ( dist < minDist )
+    {
+      minDist = dist;
+      nToReplace = n;
+    }
+  }
+
+  // remove node + replace by nToReplace
+  std::list< const SMDS_MeshNode* > nodeGroup = { nToReplace, node };
+  TListOfListOfNodes nodesToMerge( 1, nodeGroup );
+  this->MergeNodes( nodesToMerge );
+}
+
 //================================================================================
 /*!
  * \brief Create 0D elements on all nodes of the given object.
 //================================================================================
 /*!
  * \brief Create 0D elements on all nodes of the given object.
- *  \param elements - Elements on whose nodes to create 0D elements; if empty, 
+ *  \param elements - Elements on whose nodes to create 0D elements; if empty,
  *                    the all mesh is treated
  *  \param all0DElems - returns all 0D elements found or created on nodes of \a elements
  *  \param duplicateElements - to add one more 0D element to a node or not
  *                    the all mesh is treated
  *  \param all0DElems - returns all 0D elements found or created on nodes of \a elements
  *  \param duplicateElements - to add one more 0D element to a node or not
@@ -606,25 +696,6 @@ static void shiftNodesQuadTria(vector< const SMDS_MeshNode* >& aNodes)
   aNodes[5] = nd2;
 }
 
   aNodes[5] = nd2;
 }
 
-//=======================================================================
-//function : nbEdgeConnectivity
-//purpose  : return number of the edges connected with the theNode.
-//           if theEdges has connections with the other type of the
-//           elements, return -1
-//=======================================================================
-
-static int nbEdgeConnectivity(const SMDS_MeshNode* theNode)
-{
-  // SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator();
-  // int nb=0;
-  // while(elemIt->more()) {
-  //   elemIt->next();
-  //   nb++;
-  // }
-  // return nb;
-  return theNode->NbInverseElements();
-}
-
 //=======================================================================
 //function : getNodesFromTwoTria
 //purpose  : 
 //=======================================================================
 //function : getNodesFromTwoTria
 //purpose  : 
@@ -953,24 +1024,24 @@ bool getQuadrangleNodes(const SMDS_MeshNode *    theQuadNodes [],
 {
   if( tr1->NbNodes() != tr2->NbNodes() )
     return false;
 {
   if( tr1->NbNodes() != tr2->NbNodes() )
     return false;
+
   // find the 4-th node to insert into tr1
   const SMDS_MeshNode* n4 = 0;
   SMDS_ElemIteratorPtr it = tr2->nodesIterator();
   // find the 4-th node to insert into tr1
   const SMDS_MeshNode* n4 = 0;
   SMDS_ElemIteratorPtr it = tr2->nodesIterator();
-  int i=0;
-  while ( !n4 && i<3 ) {
+  for ( int i = 0; !n4 && i < 3; ++i )
+  {
     const SMDS_MeshNode * n = cast2Node( it->next() );
     const SMDS_MeshNode * n = cast2Node( it->next() );
-    i++;
     bool isDiag = ( n == theNode1 || n == theNode2 );
     if ( !isDiag )
       n4 = n;
   }
     bool isDiag = ( n == theNode1 || n == theNode2 );
     if ( !isDiag )
       n4 = n;
   }
+
   // Make an array of nodes to be in a quadrangle
   int iNode = 0, iFirstDiag = -1;
   it = tr1->nodesIterator();
   // Make an array of nodes to be in a quadrangle
   int iNode = 0, iFirstDiag = -1;
   it = tr1->nodesIterator();
-  i=0;
-  while ( i<3 ) {
+  for ( int i = 0; i < 3; ++i )
+  {
     const SMDS_MeshNode * n = cast2Node( it->next() );
     const SMDS_MeshNode * n = cast2Node( it->next() );
-    i++;
     bool isDiag = ( n == theNode1 || n == theNode2 );
     if ( isDiag ) {
       if ( iFirstDiag < 0 )
     bool isDiag = ( n == theNode1 || n == theNode2 );
     if ( isDiag ) {
       if ( iFirstDiag < 0 )
@@ -1085,6 +1156,210 @@ bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
   return true;
 }
 
   return true;
 }
 
+//=======================================================================
+//function : SplitEdge
+//purpose  : Replace each triangle bound by theNode1-theNode2 segment with
+//           two triangles by connecting a node made on the link with a node opposite to the link.
+//=======================================================================
+
+void SMESH_MeshEditor::SplitEdge (const SMDS_MeshNode * theNode1,
+                                  const SMDS_MeshNode * theNode2,
+                                  double                thePosition)
+{
+  ClearLastCreated();
+
+  SMESHDS_Mesh * mesh = GetMeshDS();
+
+  // Get triangles and segments to divide
+
+  std::vector<const SMDS_MeshNode *> diagNodes = { theNode1, theNode2 };
+  std::vector<const SMDS_MeshElement *> foundElems;
+  if ( !mesh->GetElementsByNodes( diagNodes, foundElems ) || foundElems.empty() )
+    throw SALOME_Exception( SMESH_Comment("No triangle is bound by the edge ")
+                            << theNode1->GetID() << " - " << theNode2->GetID());
+
+  SMESH_MesherHelper helper( *GetMesh() );
+
+  for ( const SMDS_MeshElement * elem : foundElems )
+  {
+    SMDSAbs_ElementType type = elem->GetType();
+    switch ( type ) {
+    case SMDSAbs_Volume:
+      throw SALOME_Exception( "Can't split an edge of a volume");
+      break;
+
+    case SMDSAbs_Face:
+      if ( elem->GetGeomType() != SMDSGeom_TRIANGLE )
+        throw SALOME_Exception( "Can't split an edge of a face of type other than triangle");
+      if ( elem->IsQuadratic() )
+      {
+        helper.SetIsQuadratic( true );
+        helper.AddTLinks( static_cast< const SMDS_MeshFace*>( elem ));
+        helper.SetIsBiQuadratic( elem->GetEntityType() == SMDSEntity_BiQuad_Triangle );
+      }
+      break;
+
+    case SMDSAbs_Edge:
+      if ( elem->IsQuadratic() )
+      {
+        helper.SetIsQuadratic( true );
+        helper.AddTLinks( static_cast< const SMDS_MeshEdge*>( elem ));
+      }
+      break;
+    default:;
+    }
+  }
+
+  // Make a new node
+
+  const SMDS_MeshNode* nodeOnLink = helper.GetMediumNode( theNode1, theNode2,/*force3d=*/false );
+
+  gp_Pnt newNodeXYZ = ( SMESH_NodeXYZ( theNode1 ) * ( 1 - thePosition ) +
+                        SMESH_NodeXYZ( theNode2 ) * thePosition );
+
+  const TopoDS_Shape& S = mesh->IndexToShape( nodeOnLink->GetShapeID() );
+  if ( !S.IsNull() && S.ShapeType() == TopAbs_FACE ) // find newNodeXYZ by UV on FACE
+  {
+    Handle(ShapeAnalysis_Surface) surface = helper.GetSurface( TopoDS::Face( S ));
+    double  tol = 100 * helper.MaxTolerance( S );
+    gp_Pnt2d uv = surface->ValueOfUV( newNodeXYZ, tol );
+    if ( surface->Gap() < SMESH_NodeXYZ( theNode1 ).Distance( theNode2 ))
+    {
+      newNodeXYZ = surface->Value( uv );
+      if ( SMDS_FacePositionPtr nPos = nodeOnLink->GetPosition())
+        nPos->SetParameters( uv.X(), uv.Y() );
+    }
+  }
+  if ( !S.IsNull() && S.ShapeType() == TopAbs_EDGE ) // find newNodeXYZ by param on EDGE
+  {
+    mesh->MoveNode( nodeOnLink, newNodeXYZ.X(), newNodeXYZ.Y(), newNodeXYZ.Z() );
+    double u = Precision::Infinite(), tol = 100 * helper.MaxTolerance( S ), distXYZ[4];
+    helper.ToFixNodeParameters( true );
+    if ( helper.CheckNodeU( TopoDS::Edge( S ), nodeOnLink, u, tol, /*force3D=*/false, distXYZ ))
+      newNodeXYZ.SetCoord( distXYZ[1], distXYZ[2], distXYZ[3] );
+  }
+  mesh->MoveNode( nodeOnLink, newNodeXYZ.X(), newNodeXYZ.Y(), newNodeXYZ.Z() );
+
+  // Split triangles and segments
+
+  std::vector<const SMDS_MeshNode *> nodes( 7 );
+  for ( const SMDS_MeshElement * elem : foundElems )
+  {
+    nodes.assign( elem->begin_nodes(), elem->end_nodes() );
+    nodes.resize( elem->NbCornerNodes() + 1 );
+    nodes.back() = nodes[0];
+
+    smIdType id = elem->GetID();
+    int shapeID = elem->GetShapeID();
+
+    const SMDS_MeshNode* centralNode = nullptr;
+    if ( elem->GetEntityType() == SMDSEntity_BiQuad_Triangle )
+      centralNode = elem->GetNode( 6 );
+
+    mesh->RemoveFreeElement( elem, /*sm=*/0, /*fromGroups=*/false );
+    if ( centralNode )
+      mesh->RemoveFreeNode( centralNode, /*sm=*/0, /*fromGroups=*/true );
+
+    for ( size_t i = 1; i < nodes.size(); ++i )
+    {
+      const SMDS_MeshNode* n1 = nodes[i-1];
+      const SMDS_MeshNode* n2 = nodes[i];
+      const SMDS_MeshElement* newElem;
+      if ( nodes.size() == 4 ) //    triangle
+      {
+        bool isDiag1 = ( n1 == theNode1 || n1 == theNode2 );
+        bool isDiag2 = ( n2 == theNode1 || n2 == theNode2 );
+        if ( isDiag1 && isDiag2 )
+          continue;
+
+        newElem = helper.AddFace( n1, n2, nodeOnLink, id );
+      }
+      else //    segment
+      {
+        newElem = helper.AddEdge( n1, nodeOnLink, id );
+      }
+      myLastCreatedElems.push_back( newElem );
+      AddToSameGroups( newElem, elem, mesh );
+      if ( shapeID )
+        mesh->SetMeshElementOnShape( newElem, shapeID );
+      id = 0;
+    }
+  }
+  return;
+}
+
+//=======================================================================
+//function : SplitFace
+//purpose  : Split a face into triangles each formed by two nodes of the 
+//           face and a new node added at the given coordinates.
+//=======================================================================
+
+void SMESH_MeshEditor::SplitFace (const SMDS_MeshElement * theFace,
+                                  double                   theX,
+                                  double                   theY,
+                                  double                   theZ )
+{
+  ClearLastCreated();
+
+  if ( !theFace )
+    throw SALOME_Exception("Null face given");
+  if ( theFace->GetType() != SMDSAbs_Face )
+    throw SALOME_Exception("Not a face given");
+
+  SMESHDS_Mesh * mesh = GetMeshDS();
+
+  SMESH_MesherHelper helper( *GetMesh() );
+  if ( theFace->IsQuadratic() )
+  {
+    helper.SetIsQuadratic( true );
+    helper.AddTLinks( static_cast< const SMDS_MeshFace*>( theFace ));
+  }
+  const TopoDS_Shape& shape = mesh->IndexToShape( theFace->GetShapeID() );
+  helper.SetSubShape( shape );
+  helper.SetElementsOnShape( true );
+
+  // Make a new node
+
+  const SMDS_MeshNode* centralNode = nullptr;
+  if (      theFace->GetEntityType() == SMDSEntity_BiQuad_Triangle )
+    centralNode = theFace->GetNode( 6 );
+  else if ( theFace->GetEntityType() == SMDSEntity_BiQuad_Quadrangle )
+    centralNode = theFace->GetNode( 8 );
+
+  if ( centralNode )
+  {
+    helper.SetIsBiQuadratic( true );
+    mesh->MoveNode( centralNode, theX, theY, theZ );
+  }
+  else
+    centralNode = helper.AddNode( theX, theY, theZ );
+
+
+  // Split theFace
+
+  std::vector<const SMDS_MeshNode *> nodes( theFace->NbNodes() + 1 );
+  nodes.assign( theFace->begin_nodes(), theFace->end_nodes() );
+  nodes.resize( theFace->NbCornerNodes() + 1 );
+  nodes.back() = nodes[0];
+
+  smIdType id = theFace->GetID();
+  int shapeID = theFace->GetShapeID();
+
+  mesh->RemoveFreeElement( theFace, /*sm=*/0, /*fromGroups=*/false );
+
+  for ( size_t i = 1; i < nodes.size(); ++i )
+  {
+    const SMDS_MeshElement* newElem = helper.AddFace( nodes[i-1], nodes[i], centralNode, id );
+
+    myLastCreatedElems.push_back( newElem );
+    AddToSameGroups( newElem, theFace, mesh );
+    if ( shapeID )
+      mesh->SetMeshElementOnShape( newElem, shapeID );
+    id = 0;
+  }
+  return;
+}
+
 //=======================================================================
 //function : Reorient
 //purpose  : Reverse theElement orientation
 //=======================================================================
 //function : Reorient
 //purpose  : Reverse theElement orientation
@@ -1112,19 +1387,37 @@ bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
       MESSAGE("Warning: bad volumic element");
       return false;
     }
       MESSAGE("Warning: bad volumic element");
       return false;
     }
-    const int nbFaces = aPolyedre->NbFaces();
+    SMDS_VolumeTool vTool( aPolyedre );
+    const int nbFaces = vTool.NbFaces();
+    vector<int> quantities( nbFaces );
     vector<const SMDS_MeshNode *> poly_nodes;
     vector<const SMDS_MeshNode *> poly_nodes;
-    vector<int> quantities (nbFaces);
 
 
-    // reverse each face of the polyedre
-    for (int iface = 1; iface <= nbFaces; iface++) {
-      int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
-      quantities[iface - 1] = nbFaceNodes;
+    // check if all facets are oriented equally
+    bool sameOri = true;
+    vector<int>& facetOri = quantities; // keep orientation in quantities so far
+    for (int iface = 0; iface < nbFaces; iface++)
+    {
+      facetOri[ iface ] = vTool.IsFaceExternal( iface );
+      if ( facetOri[ iface ] != facetOri[ 0 ])
+        sameOri = false;
+    }
 
 
-      for (inode = nbFaceNodes; inode >= 1; inode--) {
-        const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
-        poly_nodes.push_back(curNode);
-      }
+    // reverse faces of the polyhedron
+    int neededOri = sameOri ? 1 - facetOri[0] : 1;
+    poly_nodes.reserve( vTool.NbNodes() );
+    for ( int iface = 0; iface < nbFaces; iface++ )
+    {
+      int             nbFaceNodes = vTool.NbFaceNodes( iface );
+      const SMDS_MeshNode** nodes = vTool.GetFaceNodes( iface );
+      bool toReverse = ( facetOri[ iface ] != neededOri );
+
+      quantities[ iface ] = nbFaceNodes;
+
+      if ( toReverse )
+        for ( int inode = nbFaceNodes - 1; inode >= 0; inode-- )
+          poly_nodes.push_back( nodes[ inode ]);
+      else
+        poly_nodes.insert( poly_nodes.end(), nodes, nodes + nbFaceNodes );
     }
     return GetMeshDS()->ChangePolyhedronNodes( theElem, poly_nodes, quantities );
   }
     }
     return GetMeshDS()->ChangePolyhedronNodes( theElem, poly_nodes, quantities );
   }
@@ -1149,69 +1442,88 @@ bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
 /*!
  * \brief Reorient faces.
  * \param theFaces - the faces to reorient. If empty the whole mesh is meant
 /*!
  * \brief Reorient faces.
  * \param theFaces - the faces to reorient. If empty the whole mesh is meant
- * \param theDirection - desired direction of normal of \a theFace
- * \param theFace - one of \a theFaces that should be oriented according to
- *        \a theDirection and whose orientation defines orientation of other faces
+ * \param theDirection - desired direction of normal of \a theRefFaces.
+ *        It can be (0,0,0) in order to keep orientation of \a theRefFaces.
+ * \param theRefFaces - correctly oriented faces whose orientation defines
+ *        orientation of other faces.
  * \return number of reoriented faces.
  */
 //================================================================================
 
  * \return number of reoriented faces.
  */
 //================================================================================
 
-int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
-                                  const gp_Dir&            theDirection,
-                                  const SMDS_MeshElement * theFace)
+int SMESH_MeshEditor::Reorient2D( TIDSortedElemSet &  theFaces,
+                                  const gp_Vec&       theDirection,
+                                  TIDSortedElemSet &  theRefFaces,
+                                  bool                theAllowNonManifold )
 {
   int nbReori = 0;
 {
   int nbReori = 0;
-  if ( !theFace || theFace->GetType() != SMDSAbs_Face ) return nbReori;
 
   if ( theFaces.empty() )
   {
 
   if ( theFaces.empty() )
   {
-    SMDS_FaceIteratorPtr fIt = GetMeshDS()->facesIterator(/*idInceasingOrder=true*/);
+    SMDS_FaceIteratorPtr fIt = GetMeshDS()->facesIterator();
     while ( fIt->more() )
       theFaces.insert( theFaces.end(), fIt->next() );
     while ( fIt->more() )
       theFaces.insert( theFaces.end(), fIt->next() );
+
+    if ( theFaces.empty() )
+      return nbReori;
   }
 
   }
 
-  // orient theFace according to theDirection
-  gp_XYZ normal;
-  SMESH_MeshAlgos::FaceNormal( theFace, normal, /*normalized=*/false );
-  if ( normal * theDirection.XYZ() < 0 )
-    nbReori += Reorient( theFace );
+  // orient theRefFaces according to theDirection
+  if ( theDirection.X() != 0 || theDirection.Y() != 0 || theDirection.Z() != 0 )
+    for ( const SMDS_MeshElement* refFace : theRefFaces )
+    {
+      gp_XYZ normal;
+      SMESH_MeshAlgos::FaceNormal( refFace, normal, /*normalized=*/false );
+      if ( normal * theDirection.XYZ() < 0 )
+        nbReori += Reorient( refFace );
+    }
 
 
-  // Orient other faces
+  // mark reference faces
+  GetMeshDS()->SetAllCellsNotMarked();
+  for ( const SMDS_MeshElement* refFace : theRefFaces )
+    refFace->setIsMarked( true );
 
 
-  set< const SMDS_MeshElement* > startFaces, visitedFaces;
-  TIDSortedElemSet avoidSet;
-  set< SMESH_TLink > checkedLinks;
-  pair< set< SMESH_TLink >::iterator, bool > linkIt_isNew;
+  // erase reference faces from theFaces
+  for ( TIDSortedElemSet::iterator fIt = theFaces.begin(); fIt != theFaces.end(); )
+    if ( (*fIt)->isMarked() )
+      fIt = theFaces.erase( fIt );
+    else
+      ++fIt;
+
+  if ( theRefFaces.empty() )
+  {
+    theRefFaces.insert( *theFaces.begin() );
+    theFaces.erase( theFaces.begin() );
+  }
 
 
-  if ( theFaces.size() > 1 )// leave 1 face to prevent finding not selected faces
-    theFaces.erase( theFace );
-  startFaces.insert( theFace );
+  // Orient theFaces
+
+  // if ( theFaces.size() > 1 )// leave 1 face to prevent finding not selected faces
+  //   theFaces.erase( theFace );
 
   int nodeInd1, nodeInd2;
 
   int nodeInd1, nodeInd2;
-  const SMDS_MeshElement*           otherFace;
+  const SMDS_MeshElement*           refFace, *otherFace;
   vector< const SMDS_MeshElement* > facesNearLink;
   vector< std::pair< int, int > >   nodeIndsOfFace;
   vector< const SMDS_MeshElement* > facesNearLink;
   vector< std::pair< int, int > >   nodeIndsOfFace;
+  TIDSortedElemSet                  avoidSet, emptySet;
+  NCollection_Map< SMESH_TLink, SMESH_TLink > checkedLinks;
 
 
-  set< const SMDS_MeshElement* >::iterator startFace = startFaces.begin();
-  while ( !startFaces.empty() )
+  while ( !theRefFaces.empty() )
   {
   {
-    startFace = startFaces.begin();
-    theFace = *startFace;
-    startFaces.erase( startFace );
-    if ( !visitedFaces.insert( theFace ).second )
-      continue;
+    auto refFaceIt = theRefFaces.begin();
+    refFace = *refFaceIt;
+    theRefFaces.erase( refFaceIt );
 
     avoidSet.clear();
 
     avoidSet.clear();
-    avoidSet.insert(theFace);
+    avoidSet.insert( refFace );
 
 
-    NLink link( theFace->GetNode( 0 ), (SMDS_MeshNode *) 0 );
+    NLink link( refFace->GetNode( 0 ), nullptr );
 
 
-    const int nbNodes = theFace->NbCornerNodes();
-    for ( int i = 0; i < nbNodes; ++i ) // loop on links of theFace
+    const int nbNodes = refFace->NbCornerNodes();
+    for ( int i = 0; i < nbNodes; ++i ) // loop on links of refFace
     {
     {
-      link.second = theFace->GetNode(( i+1 ) % nbNodes );
-      linkIt_isNew = checkedLinks.insert( link );
-      if ( !linkIt_isNew.second )
+      link.second = refFace->GetNode(( i+1 ) % nbNodes );
+      bool isLinkVisited = checkedLinks.Contains( link );
+      if ( isLinkVisited )
       {
         // link has already been checked and won't be encountered more
         // if the group (theFaces) is manifold
       {
         // link has already been checked and won't be encountered more
         // if the group (theFaces) is manifold
@@ -1219,28 +1531,41 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
       }
       else
       {
       }
       else
       {
+        checkedLinks.Add( link );
+
         facesNearLink.clear();
         nodeIndsOfFace.clear();
         facesNearLink.clear();
         nodeIndsOfFace.clear();
+        TIDSortedElemSet::iterator objFaceIt = theFaces.end();
+
         while (( otherFace = SMESH_MeshAlgos::FindFaceInSet( link.first, link.second,
         while (( otherFace = SMESH_MeshAlgos::FindFaceInSet( link.first, link.second,
-                                                             theFaces, avoidSet,
+                                                             emptySet, avoidSet,
                                                              &nodeInd1, &nodeInd2 )))
                                                              &nodeInd1, &nodeInd2 )))
-          if ( otherFace != theFace)
+        {
+          if (( otherFace->isMarked() ) || // ref face
+              (( objFaceIt = theFaces.find( otherFace )) != theFaces.end() )) // object face
           {
             facesNearLink.push_back( otherFace );
             nodeIndsOfFace.push_back( make_pair( nodeInd1, nodeInd2 ));
           {
             facesNearLink.push_back( otherFace );
             nodeIndsOfFace.push_back( make_pair( nodeInd1, nodeInd2 ));
-            avoidSet.insert( otherFace );
           }
           }
+          avoidSet.insert( otherFace );
+        }
         if ( facesNearLink.size() > 1 )
         {
           // NON-MANIFOLD mesh shell !
         if ( facesNearLink.size() > 1 )
         {
           // NON-MANIFOLD mesh shell !
-          // select a face most co-directed with theFace,
+          if ( !theAllowNonManifold )
+          {
+            throw SALOME_Exception("Non-manifold topology of groups");
+          }
+          // select a face most co-directed with refFace,
           // other faces won't be visited this time
           gp_XYZ NF, NOF;
           // other faces won't be visited this time
           gp_XYZ NF, NOF;
-          SMESH_MeshAlgos::FaceNormal( theFace, NF, /*normalized=*/false );
+          SMESH_MeshAlgos::FaceNormal( refFace, NF, /*normalized=*/false );
           double proj, maxProj = -1;
           double proj, maxProj = -1;
-          for ( size_t i = 0; i < facesNearLink.size(); ++i ) {
+          for ( size_t i = 0; i < facesNearLink.size(); ++i )
+          {
             SMESH_MeshAlgos::FaceNormal( facesNearLink[i], NOF, /*normalized=*/false );
             SMESH_MeshAlgos::FaceNormal( facesNearLink[i], NOF, /*normalized=*/false );
-            if (( proj = Abs( NF * NOF )) > maxProj ) {
+            if (( proj = Abs( NF * NOF )) > maxProj )
+            {
               maxProj = proj;
               otherFace = facesNearLink[i];
               nodeInd1  = nodeIndsOfFace[i].first;
               maxProj = proj;
               otherFace = facesNearLink[i];
               nodeInd1  = nodeIndsOfFace[i].first;
@@ -1248,9 +1573,9 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
             }
           }
           // not to visit rejected faces
             }
           }
           // not to visit rejected faces
-          for ( size_t i = 0; i < facesNearLink.size(); ++i )
-            if ( facesNearLink[i] != otherFace && theFaces.size() > 1 )
-              visitedFaces.insert( facesNearLink[i] );
+          // for ( size_t i = 0; i < facesNearLink.size(); ++i )
+          //   if ( facesNearLink[i] != otherFace && theFaces.size() > 1 )
+          //     visitedFaces.insert( facesNearLink[i] );
         }
         else if ( facesNearLink.size() == 1 )
         {
         }
         else if ( facesNearLink.size() == 1 )
         {
@@ -1258,20 +1583,36 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
           nodeInd1  = nodeIndsOfFace.back().first;
           nodeInd2  = nodeIndsOfFace.back().second;
         }
           nodeInd1  = nodeIndsOfFace.back().first;
           nodeInd2  = nodeIndsOfFace.back().second;
         }
-        if ( otherFace && otherFace != theFace)
+        if ( otherFace )
         {
         {
-          // link must be reverse in otherFace if orientation to otherFace
-          // is same as that of theFace
-          if ( abs(nodeInd2-nodeInd1) == 1 ? nodeInd2 > nodeInd1 : nodeInd1 > nodeInd2 )
+          // link must be reverse in otherFace if orientation of otherFace
+          // is same as that of refFace
+          if ( abs( nodeInd2 - nodeInd1 ) == 1 ? nodeInd2 > nodeInd1 : nodeInd1 > nodeInd2 )
           {
           {
+            if ( otherFace->isMarked() )
+              throw SALOME_Exception("Different orientation of reference faces");
             nbReori += Reorient( otherFace );
           }
             nbReori += Reorient( otherFace );
           }
-          startFaces.insert( otherFace );
+          if ( !otherFace->isMarked() )
+          {
+            theRefFaces.insert( otherFace );
+            if ( objFaceIt != theFaces.end() )
+              theFaces.erase( objFaceIt );
+          }
         }
       }
         }
       }
-      std::swap( link.first, link.second ); // reverse the link
+      link.first = link.second; // reverse the link
+
+    } // loop on links of refFace
+
+    if ( theRefFaces.empty() && !theFaces.empty() )
+    {
+      theRefFaces.insert( *theFaces.begin() );
+      theFaces.erase( theFaces.begin() );
     }
     }
-  }
+
+  } // while ( !theRefFaces.empty() )
+
   return nbReori;
 }
 
   return nbReori;
 }
 
@@ -1494,9 +1835,12 @@ void SMESH_MeshEditor::QuadTo4Tri (TIDSortedElemSet & theElems)
   SMESH_MesherHelper helper( *GetMesh() );
   helper.SetElementsOnShape( true );
 
   SMESH_MesherHelper helper( *GetMesh() );
   helper.SetElementsOnShape( true );
 
-  SMDS_ElemIteratorPtr faceIt;
-  if ( theElems.empty() ) faceIt = GetMeshDS()->elementsIterator(SMDSAbs_Face);
-  else                    faceIt = SMESHUtils::elemSetIterator( theElems );
+  // get standalone groups of faces
+  vector< SMDS_MeshGroup* > allFaceGroups, faceGroups;
+  for ( SMESHDS_GroupBase* grBase : GetMeshDS()->GetGroups() )
+    if ( SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( grBase ))
+      if ( group->GetType() == SMDSAbs_Face && !group->IsEmpty() )
+        allFaceGroups.push_back( & group->SMDSGroup() );
 
   bool   checkUV;
   gp_XY  uv [9]; uv[8] = gp_XY(0,0);
 
   bool   checkUV;
   gp_XY  uv [9]; uv[8] = gp_XY(0,0);
@@ -1507,6 +1851,10 @@ void SMESH_MeshEditor::QuadTo4Tri (TIDSortedElemSet & theElems)
   Handle(Geom_Surface)           surface;
   TopLoc_Location                loc;
 
   Handle(Geom_Surface)           surface;
   TopLoc_Location                loc;
 
+  SMDS_ElemIteratorPtr faceIt;
+  if ( theElems.empty() ) faceIt = GetMeshDS()->elementsIterator(SMDSAbs_Face);
+  else                    faceIt = SMESHUtils::elemSetIterator( theElems );
+
   while ( faceIt->more() )
   {
     const SMDS_MeshElement* quad = faceIt->next();
   while ( faceIt->more() )
   {
     const SMDS_MeshElement* quad = faceIt->next();
@@ -1583,13 +1931,19 @@ void SMESH_MeshEditor::QuadTo4Tri (TIDSortedElemSet & theElems)
       myLastCreatedNodes.push_back( nCentral );
     }
 
       myLastCreatedNodes.push_back( nCentral );
     }
 
-    // create 4 triangles
-
     helper.SetIsQuadratic  ( nodes.size() > 4 );
     helper.SetIsBiQuadratic( nodes.size() == 9 );
     if ( helper.GetIsQuadratic() )
       helper.AddTLinks( static_cast< const SMDS_MeshFace*>( quad ));
 
     helper.SetIsQuadratic  ( nodes.size() > 4 );
     helper.SetIsBiQuadratic( nodes.size() == 9 );
     if ( helper.GetIsQuadratic() )
       helper.AddTLinks( static_cast< const SMDS_MeshFace*>( quad ));
 
+    // select groups to update
+    faceGroups.clear();
+    for ( SMDS_MeshGroup* group : allFaceGroups )
+      if ( group->Remove( quad ))
+        faceGroups.push_back( group );
+
+    // create 4 triangles
+
     GetMeshDS()->RemoveFreeElement( quad, subMeshDS, /*fromGroups=*/false );
 
     for ( int i = 0; i < 4; ++i )
     GetMeshDS()->RemoveFreeElement( quad, subMeshDS, /*fromGroups=*/false );
 
     for ( int i = 0; i < 4; ++i )
@@ -1597,8 +1951,9 @@ void SMESH_MeshEditor::QuadTo4Tri (TIDSortedElemSet & theElems)
       SMDS_MeshElement* tria = helper.AddFace( nodes[ i ],
                                                nodes[(i+1)%4],
                                                nCentral );
       SMDS_MeshElement* tria = helper.AddFace( nodes[ i ],
                                                nodes[(i+1)%4],
                                                nCentral );
-      ReplaceElemInGroups( tria, quad, GetMeshDS() );
       myLastCreatedElems.push_back( tria );
       myLastCreatedElems.push_back( tria );
+      for ( SMDS_MeshGroup* group : faceGroups )
+        group->Add( tria );
     }
   }
 }
     }
   }
 }
@@ -1784,6 +2139,17 @@ namespace
     TSplitMethod( int nbTet=0, const int* conn=0, bool addNode=false)
       : _nbSplits(nbTet), _nbCorners(4), _connectivity(conn), _baryNode(addNode), _ownConn(false) {}
     ~TSplitMethod() { if ( _ownConn ) delete [] _connectivity; _connectivity = 0; }
     TSplitMethod( int nbTet=0, const int* conn=0, bool addNode=false)
       : _nbSplits(nbTet), _nbCorners(4), _connectivity(conn), _baryNode(addNode), _ownConn(false) {}
     ~TSplitMethod() { if ( _ownConn ) delete [] _connectivity; _connectivity = 0; }
+    TSplitMethod(const TSplitMethod &splitMethod)
+      : _nbSplits(splitMethod._nbSplits),
+        _nbCorners(splitMethod._nbCorners),
+        _baryNode(splitMethod._baryNode),
+        _ownConn(splitMethod._ownConn),
+        _faceBaryNode(splitMethod._faceBaryNode)
+    {
+      _connectivity = splitMethod._connectivity;
+      const_cast<TSplitMethod&>(splitMethod)._connectivity = nullptr;
+      const_cast<TSplitMethod&>(splitMethod)._ownConn = false;
+    }
     bool hasFacet( const TTriangleFacet& facet ) const
     {
       if ( _nbCorners == 4 )
     bool hasFacet( const TTriangleFacet& facet ) const
     {
       if ( _nbCorners == 4 )
@@ -1903,7 +2269,7 @@ namespace
       }
       for ( int variant = 0; variant < nbVariants && method._nbSplits == 0; ++variant )
       {
       }
       for ( int variant = 0; variant < nbVariants && method._nbSplits == 0; ++variant )
       {
-        // check method compliancy with adjacent tetras,
+        // check method compliance with adjacent tetras,
         // all found splits must be among facets of tetras described by this method
         method = TSplitMethod( nbTet, connVariants[variant] );
         if ( hasAdjacentSplits && method._nbSplits > 0 )
         // all found splits must be among facets of tetras described by this method
         method = TSplitMethod( nbTet, connVariants[variant] );
         if ( hasAdjacentSplits && method._nbSplits > 0 )
@@ -2034,6 +2400,8 @@ namespace
                                     const int        methodFlags,
                                     const int        facetToSplit)
   {
                                     const int        methodFlags,
                                     const int        facetToSplit)
   {
+    TSplitMethod method;
+
     // order of facets in HEX according to SMDS_VolumeTool::Hexa_F :
     // B, T, L, B, R, F
     const int iF = ( facetToSplit < 2 ) ? 0 : 1 + ( facetToSplit-2 ) % 2; // [0,1,2]
     // order of facets in HEX according to SMDS_VolumeTool::Hexa_F :
     // B, T, L, B, R, F
     const int iF = ( facetToSplit < 2 ) ? 0 : 1 + ( facetToSplit-2 ) % 2; // [0,1,2]
@@ -2064,12 +2432,12 @@ namespace
         to4methods[iF]._nbSplits  = 4;
         to4methods[iF]._nbCorners = 6;
       }
         to4methods[iF]._nbSplits  = 4;
         to4methods[iF]._nbCorners = 6;
       }
-      return to4methods[iF];
+      method = to4methods[iF];
+      to4methods[iF]._connectivity = method._connectivity; // as copy ctor resets _connectivity
+      return method;
     }
     // else if ( methodFlags == HEXA_TO_2_PRISMS )
 
     }
     // else if ( methodFlags == HEXA_TO_2_PRISMS )
 
-    TSplitMethod method;
-
     const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
 
     const int nbVariants = 2, nbSplits = 2;
     const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
 
     const int nbVariants = 2, nbSplits = 2;
@@ -2102,7 +2470,7 @@ namespace
       // there are adjacent prism
       for ( int variant = 0; variant < nbVariants; ++variant )
       {
       // there are adjacent prism
       for ( int variant = 0; variant < nbVariants; ++variant )
       {
-        // check method compliancy with adjacent prisms,
+        // check method compliance with adjacent prisms,
         // the found prism facets must be among facets of prisms described by current method
         method._nbSplits     = nbSplits;
         method._nbCorners    = 6;
         // the found prism facets must be among facets of prisms described by current method
         method._nbSplits     = nbSplits;
         method._nbCorners    = 6;
@@ -2181,7 +2549,7 @@ namespace
    */
   //=======================================================================
 
    */
   //=======================================================================
 
-  struct TVolumeFaceKey: pair< pair< int, int>, pair< int, int> >
+  struct TVolumeFaceKey: pair< pair< smIdType, smIdType>, pair< smIdType, smIdType> >
   {
     TVolumeFaceKey( SMDS_VolumeTool& vol, int iF )
     {
   {
     TVolumeFaceKey( SMDS_VolumeTool& vol, int iF )
     {
@@ -3095,10 +3463,10 @@ public:
     :myMesh( theMesh ), myMaxID( theMesh->MaxNodeID() + 1)
   {}
 
     :myMesh( theMesh ), myMaxID( theMesh->MaxNodeID() + 1)
   {}
 
-  long GetLinkID (const SMDS_MeshNode * n1,
+  smIdType GetLinkID (const SMDS_MeshNode * n1,
                   const SMDS_MeshNode * n2) const
   {
                   const SMDS_MeshNode * n2) const
   {
-    return ( Min(n1->GetID(),n2->GetID()) * myMaxID + Max(n1->GetID(),n2->GetID()));
+    return ( std::min(n1->GetID(),n2->GetID()) * myMaxID + std::max(n1->GetID(),n2->GetID()));
   }
 
   bool GetNodes (const long             theLinkID,
   }
 
   bool GetNodes (const long             theLinkID,
@@ -3605,7 +3973,7 @@ static bool getClosestUV (Extrema_GenExtPS& projector,
 {
   projector.Perform( point );
   if ( projector.IsDone() ) {
 {
   projector.Perform( point );
   if ( projector.IsDone() ) {
-    double u, v, minVal = DBL_MAX;
+    double u = 0, v = 0, minVal = DBL_MAX;
     for ( int i = projector.NbExt(); i > 0; i-- )
       if ( projector.SquareDistance( i ) < minVal ) {
         minVal = projector.SquareDistance( i );
     for ( int i = projector.NbExt(); i > 0; i-- )
       if ( projector.SquareDistance( i ) < minVal ) {
         minVal = projector.SquareDistance( i );
@@ -5183,6 +5551,7 @@ SMESH_MeshEditor::RotationSweep(TIDSortedElemSet   theElemSets[2],
 SMESH_MeshEditor::ExtrusParam::ExtrusParam( const gp_Vec&            theStep,
                                             const int                theNbSteps,
                                             const std::list<double>& theScales,
 SMESH_MeshEditor::ExtrusParam::ExtrusParam( const gp_Vec&            theStep,
                                             const int                theNbSteps,
                                             const std::list<double>& theScales,
+                                            const std::list<double>& theAngles,
                                             const gp_XYZ*            theBasePoint,
                                             const int                theFlags,
                                             const double             theTolerance):
                                             const gp_XYZ*            theBasePoint,
                                             const int                theFlags,
                                             const double             theTolerance):
@@ -5197,32 +5566,52 @@ SMESH_MeshEditor::ExtrusParam::ExtrusParam( const gp_Vec&            theStep,
   for (int i=1; i<=theNbSteps; i++ )
     mySteps->Append( stepSize );
 
   for (int i=1; i<=theNbSteps; i++ )
     mySteps->Append( stepSize );
 
-  int nbScales = theScales.size();
-  if ( nbScales > 0 )
+  if ( !theScales.empty() )
   {
   {
-    if ( IsLinearVariation() && nbScales < theNbSteps )
+    if ( IsScaleVariation() && (int)theScales.size() < theNbSteps )
+      linearScaleVariation( theNbSteps, const_cast< std::list<double>& >( theScales ));
+
+    // add medium scales
+    std::list<double>::const_iterator s2 = theScales.begin(), s1 = s2++;
+    myScales.reserve( theNbSteps * 2 );
+    myScales.push_back( 0.5 * ( *s1 + 1. ));
+    myScales.push_back( *s1 );
+    for ( ; s2 != theScales.end(); s1 = s2++ )
     {
     {
-      myScales.reserve( theNbSteps );
-      std::list<double>::const_iterator scale = theScales.begin();
-      double prevScale = 1.0;
-      for ( int iSc = 1; scale != theScales.end(); ++scale, ++iSc )
-      {
-        int      iStep = int( iSc / double( nbScales ) * theNbSteps + 0.5 );
-        int    stDelta = Max( 1, iStep - myScales.size());
-        double scDelta = ( *scale - prevScale ) / stDelta;
-        for ( int iStep = 0; iStep < stDelta; ++iStep )
-        {
-          myScales.push_back( prevScale + scDelta );
-          prevScale = myScales.back();
-        }
-        prevScale = *scale;
-      }
+      myScales.push_back( 0.5 * ( *s1 + *s2 ));
+      myScales.push_back( *s2 );
     }
     }
-    else
+  }
+
+  if ( !theAngles.empty() )
+  {
+    std::list<double>& angles = const_cast< std::list<double>& >( theAngles );
+    if ( IsAngleVariation() && (int)theAngles.size() < theNbSteps )
+      linearAngleVariation( theNbSteps, angles );
+
+    // accumulate angles
+    double angle = 0;
+    int nbAngles = 0;
+    std::list<double>::iterator a1 = angles.begin(), a2;
+    for ( ; a1 != angles.end(); ++a1, ++nbAngles )
+    {
+      angle += *a1;
+      *a1 = angle;
+    }
+    while ( nbAngles++ < theNbSteps )
+      angles.push_back( angles.back() );
+
+    // add medium angles
+    a2 = angles.begin(), a1 = a2++;
+    myAngles.push_back( 0.5 * *a1 );
+    myAngles.push_back( *a1 );
+    for ( ; a2 != angles.end(); a1 = a2++ )
     {
     {
-      myScales.assign( theScales.begin(), theScales.end() );
+      myAngles.push_back( 0.5 * ( *a1 + *a2 ));
+      myAngles.push_back( *a2 );
     }
   }
     }
   }
+
   if ( theBasePoint )
   {
     myBaseP = *theBasePoint;
   if ( theBasePoint )
   {
     myBaseP = *theBasePoint;
@@ -5293,6 +5682,41 @@ SMESH_MeshEditor::ExtrusParam::ExtrusParam( const double theStepSize,
   }
 }
 
   }
 }
 
+//=======================================================================
+//function : ExtrusParam
+//purpose  : for extrusion along path
+//=======================================================================
+
+SMESH_MeshEditor::ExtrusParam::ExtrusParam( const std::vector< PathPoint >& thePoints,
+                                            const gp_Pnt*                   theBasePoint,
+                                            const std::list<double>&        theScales,
+                                            const bool                      theMakeGroups )
+  : myBaseP( Precision::Infinite(), 0, 0 ),
+    myFlags( EXTRUSION_FLAG_BOUNDARY | ( theMakeGroups ? EXTRUSION_FLAG_GROUPS : 0 )),
+    myPathPoints( thePoints )
+{
+  if ( theBasePoint )
+  {
+    myBaseP = theBasePoint->XYZ();
+  }
+
+  if ( !theScales.empty() )
+  {
+    // add medium scales
+    std::list<double>::const_iterator s2 = theScales.begin(), s1 = s2++;
+    myScales.reserve( thePoints.size() * 2 );
+    myScales.push_back( 0.5 * ( 1. + *s1 ));
+    myScales.push_back( *s1 );
+    for ( ; s2 != theScales.end(); s1 = s2++ )
+    {
+      myScales.push_back( 0.5 * ( *s1 + *s2 ));
+      myScales.push_back( *s2 );
+    }
+  }
+
+  myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesAlongTrack;
+}
+
 //=======================================================================
 //function : ExtrusParam::SetElementsToUse
 //purpose  : stores elements to use for extrusion by normal, depending on
 //=======================================================================
 //function : ExtrusParam::SetElementsToUse
 //purpose  : stores elements to use for extrusion by normal, depending on
@@ -5398,38 +5822,37 @@ makeNodesByDir( SMESHDS_Mesh*                     mesh,
     newNodes.push_back( newNode );
   }
 
     newNodes.push_back( newNode );
   }
 
-  if ( !myScales.empty() )
+  if ( !myScales.empty() || !myAngles.empty() )
   {
   {
-    if ( makeMediumNodes && myMediumScales.empty() )
-    {
-      myMediumScales.resize( myScales.size() );
-      double prevFactor = 1.;
-      for ( size_t i = 0; i < myScales.size(); ++i )
-      {
-        myMediumScales[i] = 0.5 * ( prevFactor + myScales[i] );
-        prevFactor = myScales[i];
-      }
-    }
-    typedef std::vector<double>::iterator ScaleIt;
-    ScaleIt scales[] = { myScales.begin(), myMediumScales.begin() };
+    gp_XYZ  center = myBaseP;
+    gp_Ax1  ratationAxis( center, myDir );
+    gp_Trsf rotation;
 
 
-    size_t iSc = 0, nbScales = myScales.size() + myMediumScales.size();
-
-    gp_XYZ center = myBaseP;
     std::list<const SMDS_MeshNode*>::iterator nIt = newNodes.begin();
     std::list<const SMDS_MeshNode*>::iterator nIt = newNodes.begin();
-    size_t iN  = 0;
-    for ( beginStepIter( makeMediumNodes ); moreSteps() && ( iN < nbScales ); ++nIt, ++iN )
+    size_t i = !makeMediumNodes;
+    for ( beginStepIter( makeMediumNodes );
+          moreSteps();
+          ++nIt, i += 1 + !makeMediumNodes )
     {
       center += myDir.XYZ() * nextStep();
 
     {
       center += myDir.XYZ() * nextStep();
 
-      iSc += int( makeMediumNodes );
-      ScaleIt& scale = scales[ iSc % 2 ];
-      
       gp_XYZ xyz = SMESH_NodeXYZ( *nIt );
       gp_XYZ xyz = SMESH_NodeXYZ( *nIt );
-      xyz = ( *scale * ( xyz - center )) + center;
-      mesh->MoveNode( *nIt, xyz.X(), xyz.Y(), xyz.Z() );
-
-      ++scale;
+      bool moved = false;
+      if ( i < myScales.size() )
+      {
+        xyz = ( myScales[i] * ( xyz - center )) + center;
+        moved = true;
+      }
+      if ( !myAngles.empty() )
+      {
+        rotation.SetRotation( ratationAxis, myAngles[i] );
+        rotation.Transforms( xyz );
+        moved = true;
+      }
+      if ( moved )
+        mesh->MoveNode( *nIt, xyz.X(), xyz.Y(), xyz.Z() );
+      else
+        break;
     }
   }
   return nbNodes;
     }
   }
   return nbNodes;
@@ -5587,13 +6010,107 @@ makeNodesByNormal2D( SMESHDS_Mesh*                     mesh,
 //=======================================================================
 
 int SMESH_MeshEditor::ExtrusParam::
 //=======================================================================
 
 int SMESH_MeshEditor::ExtrusParam::
-makeNodesByNormal1D( SMESHDS_Mesh*                     mesh,
+makeNodesByNormal1D( SMESHDS_Mesh*                     /*mesh*/,
+                     const SMDS_MeshNode*              /*srcNode*/,
+                     std::list<const SMDS_MeshNode*> & /*newNodes*/,
+                     const bool                        /*makeMediumNodes*/)
+{
+  throw SALOME_Exception("Extrusion 1D by Normal not implemented");
+  return 0;
+}
+
+//=======================================================================
+//function : ExtrusParam::makeNodesAlongTrack
+//purpose  : create nodes for extrusion along path
+//=======================================================================
+
+int SMESH_MeshEditor::ExtrusParam::
+makeNodesAlongTrack( SMESHDS_Mesh*                     mesh,
                      const SMDS_MeshNode*              srcNode,
                      std::list<const SMDS_MeshNode*> & newNodes,
                      const bool                        makeMediumNodes)
 {
                      const SMDS_MeshNode*              srcNode,
                      std::list<const SMDS_MeshNode*> & newNodes,
                      const bool                        makeMediumNodes)
 {
-  throw SALOME_Exception("Extrusion 1D by Normal not implemented");
-  return 0;
+  const Standard_Real aTolAng=1.e-4;
+
+  gp_Pnt aV0x = myBaseP;
+  gp_Pnt aPN0 = SMESH_NodeXYZ( srcNode );
+
+  const PathPoint& aPP0 = myPathPoints[0];
+  gp_Pnt aP0x = aPP0.myPnt;
+  gp_Dir aDT0x= aPP0.myTgt;
+
+  std::vector< gp_Pnt > centers;
+  centers.reserve( NbSteps() * 2 );
+
+  gp_Trsf aTrsf, aTrsfRot, aTrsfRotT1T0;
+
+  for ( size_t j = 1; j < myPathPoints.size(); ++j )
+  {
+    const PathPoint&  aPP  = myPathPoints[j];
+    const gp_Pnt&     aP1x = aPP.myPnt;
+    const gp_Dir&    aDT1x = aPP.myTgt;
+
+    // Translation
+    gp_Vec aV01x( aP0x, aP1x );
+    aTrsf.SetTranslation( aV01x );
+    gp_Pnt aV1x = aV0x.Transformed( aTrsf );
+    gp_Pnt aPN1 = aPN0.Transformed( aTrsf );
+
+    // rotation 1 [ T1,T0 ]
+    Standard_Real aAngleT1T0 = -aDT1x.Angle( aDT0x );
+    if ( fabs( aAngleT1T0 ) > aTolAng )
+    {
+      gp_Dir aDT1T0 = aDT1x ^ aDT0x;
+      aTrsfRotT1T0.SetRotation( gp_Ax1( aV1x, aDT1T0 ), aAngleT1T0 );
+
+      aPN1 = aPN1.Transformed( aTrsfRotT1T0 );
+    }
+
+    // rotation 2
+    if ( aPP.myAngle != 0. )
+    {
+      aTrsfRot.SetRotation( gp_Ax1( aV1x, aDT1x ), aPP.myAngle );
+      aPN1 = aPN1.Transformed( aTrsfRot );
+    }
+
+    // make new node
+    if ( makeMediumNodes )
+    {
+      // create additional node
+      gp_XYZ midP = 0.5 * ( aPN1.XYZ() + aPN0.XYZ() );
+      const SMDS_MeshNode* newNode = mesh->AddNode( midP.X(), midP.Y(), midP.Z() );
+      newNodes.push_back( newNode );
+
+    }
+    const SMDS_MeshNode* newNode = mesh->AddNode( aPN1.X(), aPN1.Y(), aPN1.Z() );
+    newNodes.push_back( newNode );
+
+    centers.push_back( 0.5 * ( aV0x.XYZ() + aV1x.XYZ() ));
+    centers.push_back( aV1x );
+
+    aPN0 = aPN1;
+    aP0x = aP1x;
+    aV0x = aV1x;
+    aDT0x = aDT1x;
+  }
+
+  // scale
+  if ( !myScales.empty() )
+  {
+    gp_Trsf aTrsfScale;
+    std::list<const SMDS_MeshNode*>::iterator node = newNodes.begin();
+    for ( size_t i = !makeMediumNodes;
+          i < myScales.size() && node != newNodes.end();
+          i += ( 1 + !makeMediumNodes ), ++node )
+    {
+      aTrsfScale.SetScale( centers[ i ], myScales[ i ] );
+      gp_Pnt aN = SMESH_NodeXYZ( *node );
+      gp_Pnt aP = aN.Transformed( aTrsfScale );
+      mesh->MoveNode( *node, aP.X(), aP.Y(), aP.Z() );
+    }
+  }
+
+  return myPathPoints.size() + makeMediumNodes * ( myPathPoints.size() - 2 );
 }
 
 //=======================================================================
 }
 
 //=======================================================================
@@ -5609,10 +6126,31 @@ SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet     theElems[2],
                                   const int            theFlags,
                                   const double         theTolerance)
 {
                                   const int            theFlags,
                                   const double         theTolerance)
 {
-  ExtrusParam aParams( theStep, theNbSteps, std::list<double>(), 0, theFlags, theTolerance );
+  std::list<double> dummy;
+  ExtrusParam aParams( theStep, theNbSteps, dummy, dummy, 0,
+                       theFlags, theTolerance );
   return ExtrusionSweep( theElems, aParams, newElemsMap );
 }
 
   return ExtrusionSweep( theElems, aParams, newElemsMap );
 }
 
+namespace
+{
+
+//=======================================================================
+//function : getOriFactor
+//purpose  : Return -1 or 1 depending on if order of given nodes corresponds to
+//           edge curve orientation
+//=======================================================================
+
+  double getOriFactor( const TopoDS_Edge&   edge,
+                       const SMDS_MeshNode* n1,
+                       const SMDS_MeshNode* n2,
+                       SMESH_MesherHelper&  helper)
+  {
+    double u1 = helper.GetNodeU( edge, n1, n2 );
+    double u2 = helper.GetNodeU( edge, n2, n1 );
+    return u1 < u2 ? 1. : -1.;
+  }
+}
 
 //=======================================================================
 //function : ExtrusionSweep
 
 //=======================================================================
 //function : ExtrusionSweep
@@ -5661,11 +6199,11 @@ SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet     theElemSets[2],
       newNodesItVec.reserve( nbNodes );
 
       // loop on elem nodes
       newNodesItVec.reserve( nbNodes );
 
       // loop on elem nodes
-      SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+      SMDS_NodeIteratorPtr itN = elem->nodeIterator();
       while ( itN->more() )
       {
         // check if a node has been already sweeped
       while ( itN->more() )
       {
         // check if a node has been already sweeped
-        const SMDS_MeshNode* node = cast2Node( itN->next() );
+        const SMDS_MeshNode* node = itN->next();
         TNodeOfNodeListMap::iterator nIt =
           mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
         list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
         TNodeOfNodeListMap::iterator nIt =
           mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
         list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
@@ -5699,6 +6237,11 @@ SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet     theElemSets[2],
           }
           else
           {
           }
           else
           {
+            if ( theParams.ToMakeBoundary() )
+            {
+              GetMeshDS()->Modified();
+              throw SALOME_Exception( SMESH_Comment("Can't extrude node #") << node->GetID() );
+            }
             break; // newNodesItVec will be shorter than nbNodes
           }
         }
             break; // newNodesItVec will be shorter than nbNodes
           }
         }
@@ -5726,726 +6269,184 @@ SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet     theElemSets[2],
 //=======================================================================
 SMESH_MeshEditor::Extrusion_Error
 SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet     theElements[2],
 //=======================================================================
 SMESH_MeshEditor::Extrusion_Error
 SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet     theElements[2],
-                                       SMESH_subMesh*       theTrack,
+                                       SMESH_Mesh*          theTrackMesh,
+                                       SMDS_ElemIteratorPtr theTrackIterator,
                                        const SMDS_MeshNode* theN1,
                                        const SMDS_MeshNode* theN1,
-                                       const bool           theHasAngles,
-                                       list<double>&        theAngles,
-                                       const bool           theLinearVariation,
-                                       const bool           theHasRefPoint,
-                                       const gp_Pnt&        theRefPoint,
+                                       std::list<double>&   theAngles,
+                                       const bool           theAngleVariation,
+                                       std::list<double>&   theScales,
+                                       const bool           theScaleVariation,
+                                       const gp_Pnt*        theRefPoint,
                                        const bool           theMakeGroups)
 {
   ClearLastCreated();
 
                                        const bool           theMakeGroups)
 {
   ClearLastCreated();
 
-  int aNbE;
-  std::list<double> aPrms;
-  TIDSortedElemSet::iterator itElem;
-
-  gp_XYZ aGC;
-  TopoDS_Edge aTrackEdge;
-  TopoDS_Vertex aV1, aV2;
-
-  SMDS_ElemIteratorPtr aItE;
-  SMDS_NodeIteratorPtr aItN;
-  SMDSAbs_ElementType aTypeE;
-
-  TNodeOfNodeListMap mapNewNodes;
-
   // 1. Check data
   // 1. Check data
-  aNbE = theElements[0].size() + theElements[1].size();
-  // nothing to do
-  if ( !aNbE )
+  if ( theElements[0].empty() && theElements[1].empty() )
     return EXTR_NO_ELEMENTS;
 
     return EXTR_NO_ELEMENTS;
 
-  // 1.1 Track Pattern
-  ASSERT( theTrack );
-
-  SMESHDS_SubMesh* pSubMeshDS = theTrack->GetSubMeshDS();
-  if ( !pSubMeshDS )
-    return ExtrusionAlongTrack( theElements, theTrack->GetFather(), theN1,
-                                theHasAngles, theAngles, theLinearVariation,
-                                theHasRefPoint, theRefPoint, theMakeGroups );
-
-  aItE = pSubMeshDS->GetElements();
-  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;
-  }
-
-  list<SMESH_MeshEditor_PathPoint> fullList;
-
-  const TopoDS_Shape& aS = theTrack->GetSubShape();
-  // 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 ( SMESH_Algo::isDegenerated( 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();
-      SMDS_EdgePositionPtr pEPos = pNode->GetPosition();
-      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();
-          SMDS_EdgePositionPtr pEPos = pNode->GetPosition();
-          double aT = pEPos->GetUParameter();
-          aPrms.push_back( aT );
-        }
-        list<SMESH_MeshEditor_PathPoint> LPP;
-        //Extrusion_Error err =
-        makeEdgePathPoints(aPrms, aTrackEdge,(aN1->GetID()==startNid), LPP);
-        LLPPs.push_back(LPP);
-        UsedNums.Add(k);
-        // update startN for search following edge
-        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;
-  }
-
-  return makeExtrElements(theElements, fullList, theHasAngles, theAngles, theLinearVariation,
-                          theHasRefPoint, theRefPoint, theMakeGroups);
-}
-
-
-//=======================================================================
-//function : ExtrusionAlongTrack
-//purpose  :
-//=======================================================================
-SMESH_MeshEditor::Extrusion_Error
-SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet     theElements[2],
-                                       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)
-{
-  ClearLastCreated();
-
-  int aNbE;
-  std::list<double> aPrms;
-  TIDSortedElemSet::iterator itElem;
-
-  gp_XYZ aGC;
-  TopoDS_Edge aTrackEdge;
-  TopoDS_Vertex aV1, aV2;
-
-  SMDS_ElemIteratorPtr aItE;
-  SMDS_NodeIteratorPtr aItN;
-  SMDSAbs_ElementType aTypeE;
+  ASSERT( theTrackMesh );
+  if ( ! theTrackIterator || !theTrackIterator->more() )
+    return EXTR_NO_ELEMENTS;
 
 
-  TNodeOfNodeListMap mapNewNodes;
+  // 2. Get ordered nodes
+  SMESH_MeshAlgos::TElemGroupVector branchEdges;
+  SMESH_MeshAlgos::TNodeGroupVector branchNods;
+  SMESH_MeshAlgos::Get1DBranches( theTrackIterator, branchEdges, branchNods, theN1 );
+  if ( branchEdges.empty() )
+    return EXTR_PATH_NOT_EDGE;
 
 
-  // 1. Check data
-  aNbE = theElements[0].size() + theElements[1].size();
-  // nothing to do
-  if ( !aNbE )
-    return EXTR_NO_ELEMENTS;
+  if ( branchEdges.size() > 1 )
+    return EXTR_BAD_PATH_SHAPE;
 
 
-  // 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;
-  }
-
-  list<SMESH_MeshEditor_PathPoint> fullList;
-
-  const TopoDS_Shape& aS = theTrack->GetShapeToMesh();
-
-  if ( !theTrack->HasShapeToMesh() ) {
-    //Mesh without shape
-    const SMDS_MeshNode* currentNode = NULL;
-    const SMDS_MeshNode* prevNode = theN1;
-    std::vector<const SMDS_MeshNode*> aNodesList;
-    aNodesList.push_back(theN1);
-    int nbEdges = 0, conn=0;
-    const SMDS_MeshElement* prevElem = NULL;
-    const SMDS_MeshElement* currentElem = NULL;
-    int totalNbEdges = theTrack->NbEdges();
-    SMDS_ElemIteratorPtr nIt;
-
-    //check start node
-    if( !theTrack->GetMeshDS()->Contains( theN1 )) {
-      return EXTR_BAD_STARTING_NODE;
-    }
-
-    conn = nbEdgeConnectivity(theN1);
-    if( conn != 1 )
-      return EXTR_PATH_NOT_EDGE;
-
-    aItE = theN1->GetInverseElementIterator();
-    prevElem = aItE->next();
-    currentElem = prevElem;
-    //Get all nodes
-    if(totalNbEdges == 1 ) {
-      nIt = currentElem->nodesIterator();
-      currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
-      if(currentNode == prevNode)
-        currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
-      aNodesList.push_back(currentNode);
-    } else {
-      nIt = currentElem->nodesIterator();
-      while( nIt->more() ) {
-        currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
-        if(currentNode == prevNode)
-          currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
-        aNodesList.push_back(currentNode);
-
-        //case of the closed mesh
-        if(currentNode == theN1) {
-          nbEdges++;
-          break;
-        }
+  std::vector< const SMDS_MeshNode* >&    pathNodes = branchNods[0];
+  std::vector< const SMDS_MeshElement* >& pathEdges = branchEdges[0];
+  if ( pathNodes[0] != theN1 && pathNodes[1] != theN1 )
+    return EXTR_BAD_STARTING_NODE;
 
 
-        conn = nbEdgeConnectivity(currentNode);
-        if(conn > 2) {
-          return EXTR_PATH_NOT_EDGE;
-        }else if( conn == 1 && nbEdges > 0 ) {
-          //End of the path
-          nbEdges++;
-          break;
-        }else {
-          prevNode = currentNode;
-          aItE = currentNode->GetInverseElementIterator();
-          currentElem = aItE->next();
-          if( currentElem  == prevElem)
-            currentElem = aItE->next();
-          nIt = currentElem->nodesIterator();
-          prevElem = currentElem;
-          nbEdges++;
-        }
-      }
-    }
-
-    if(nbEdges != totalNbEdges)
-      return EXTR_PATH_NOT_EDGE;
-
-    TopTools_SequenceOfShape Edges;
-    list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
-    int startNid = theN1->GetID();
-    for ( size_t i = 1; i < aNodesList.size(); i++ )
-    {
-      gp_Pnt     p1 = SMESH_NodeXYZ( aNodesList[i-1] );
-      gp_Pnt     p2 = SMESH_NodeXYZ( aNodesList[i] );
-      TopoDS_Edge e = BRepBuilderAPI_MakeEdge( p1, p2 );
-      list<SMESH_MeshEditor_PathPoint> LPP;
-      aPrms.clear();
-      makeEdgePathPoints(aPrms, e, (aNodesList[i-1]->GetID()==startNid), LPP);
-      LLPPs.push_back(LPP);
-      if ( aNodesList[i-1]->GetID() == startNid ) startNid = aNodesList[i  ]->GetID();
-      else                                        startNid = aNodesList[i-1]->GetID();
-    }
-
-    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();
-    SMESH_MeshEditor_PathPoint PP2;
-    fullList.pop_back();
-    itLLPP++;
-    for(; itLLPP!=LLPPs.end(); itLLPP++) {
-      list<SMESH_MeshEditor_PathPoint> currList = *itLLPP;
-      itPP = currList.begin();
-      PP2 = currList.front();
-      gp_Dir D1 = PP1.Tangent();
-      gp_Dir D2 = PP2.Tangent();
-      gp_Dir Dnew( 0.5 * ( D1.XYZ() + D2.XYZ() ));
-      PP1.SetTangent(Dnew);
-      fullList.push_back(PP1);
-      itPP++;
-      for(; itPP!=currList.end(); itPP++) {
-        fullList.push_back( *itPP );
-      }
-      PP1 = fullList.back();
-      fullList.pop_back();
-    }
-    fullList.push_back(PP1);
-
-  } // Sub-shape for the Pattern must be an Edge or Wire
-  else if ( aS.ShapeType() == TopAbs_EDGE )
-  {
-    aTrackEdge = TopoDS::Edge( aS );
-    // the Edge must not be degenerated
-    if ( SMESH_Algo::isDegenerated( aTrackEdge ) )
-      return EXTR_BAD_PATH_SHAPE;
-    TopExp::Vertices( aTrackEdge, aV1, aV2 );
-    const SMDS_MeshNode* aN1 = SMESH_Algo::VertexNode( aV1, pMeshDS );
-    const SMDS_MeshNode* aN2 = SMESH_Algo::VertexNode( aV2, pMeshDS );
-    // 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;
-      SMDS_EdgePositionPtr pEPos = pNode->GetPosition();
-      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( SMESH_Algo::isDegenerated(E) ) continue;
-      SMESH_subMesh* SM = theTrack->GetSubMesh(E);
-      if(SM) {
-        LSM.push_back(SM);
-        Edges.Append(E);
-      }
-    }
-    list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
-    TopoDS_Vertex aVprev;
-    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 );
-        bool aN1isOK = false, aN2isOK = false;
-        if ( aVprev.IsNull() ) {
-          // if previous vertex is not yet defined, it means that we in the beginning of wire
-          // and we have to find initial vertex corresponding to starting node theN1
-          const SMDS_MeshNode* aN1 = SMESH_Algo::VertexNode( aV1, pMeshDS );
-          const SMDS_MeshNode* aN2 = SMESH_Algo::VertexNode( aV2, pMeshDS );
-          // starting node must be aN1 or aN2
-          aN1isOK = ( aN1 && aN1 == theN1 );
-          aN2isOK = ( aN2 && aN2 == theN1 );
-        }
-        else {
-          // we have specified ending vertex of the previous edge on the previous iteration
-          // and we have just to check that it corresponds to any vertex in current segment
-          aN1isOK = aVprev.IsSame( aV1 );
-          aN2isOK = aVprev.IsSame( aV2 );
-        }
-        if ( !aN1isOK && !aN2isOK ) continue;
-        // 2. Collect parameters on the track edge
-        aPrms.clear();
-        aItN = locMeshDS->GetNodes();
-        while ( aItN->more() ) {
-          const SMDS_MeshNode*  pNode = aItN->next();
-          SMDS_EdgePositionPtr pEPos = pNode->GetPosition();
-          double aT = pEPos->GetUParameter();
-          aPrms.push_back( aT );
-        }
-        list<SMESH_MeshEditor_PathPoint> LPP;
-        //Extrusion_Error err =
-        makeEdgePathPoints(aPrms, aTrackEdge, aN1isOK, LPP);
-        LLPPs.push_back(LPP);
-        UsedNums.Add(k);
-        // update startN for search following edge
-        if ( aN1isOK ) aVprev = aV2;
-        else           aVprev = aV1;
-        break;
+  if ( theTrackMesh->NbEdges( ORDER_QUADRATIC ) > 0 )
+  {
+    // add medium nodes to pathNodes
+    std::vector< const SMDS_MeshNode* >    pathNodes2;
+    std::vector< const SMDS_MeshElement* > pathEdges2;
+    pathNodes2.reserve( pathNodes.size() * 2 );
+    pathEdges2.reserve( pathEdges.size() * 2 );
+    for ( size_t i = 0; i < pathEdges.size(); ++i )
+    {
+      pathNodes2.push_back( pathNodes[i] );
+      pathEdges2.push_back( pathEdges[i] );
+      if ( pathEdges[i]->IsQuadratic() )
+      {
+        pathNodes2.push_back( pathEdges[i]->GetNode(2) );
+        pathEdges2.push_back( pathEdges[i] );
       }
     }
       }
     }
-    list< list<SMESH_MeshEditor_PathPoint> >::iterator itLLPP = LLPPs.begin();
-    list<SMESH_MeshEditor_PathPoint>& firstList = *itLLPP;
-    fullList.splice( fullList.end(), firstList );
-
-    SMESH_MeshEditor_PathPoint PP1 = fullList.back();
-    fullList.pop_back();
-    itLLPP++;
-    for(; itLLPP!=LLPPs.end(); itLLPP++) {
-      list<SMESH_MeshEditor_PathPoint>& currList = *itLLPP;
-      SMESH_MeshEditor_PathPoint PP2 = currList.front();
-      gp_Dir D1 = PP1.Tangent();
-      gp_Dir D2 = PP2.Tangent();
-      gp_Dir Dnew( D1.XYZ() + D2.XYZ() );
-      PP1.SetTangent(Dnew);
-      fullList.push_back(PP1);
-      fullList.splice( fullList.end(), currList, ++currList.begin(), currList.end() );
-      PP1 = fullList.back();
-      fullList.pop_back();
-    }
-    // if wire not closed
-    fullList.push_back(PP1);
-    // else ???
-  }
-  else {
-    return EXTR_BAD_PATH_SHAPE;
+    pathNodes2.push_back( pathNodes.back() );
+    pathEdges.swap( pathEdges2 );
+    pathNodes.swap( pathNodes2 );
   }
 
   }
 
-  return makeExtrElements(theElements, fullList, theHasAngles, theAngles, theLinearVariation,
-                          theHasRefPoint, theRefPoint, theMakeGroups);
-}
+  // 3. Get path data at pathNodes
 
 
+  std::vector< ExtrusParam::PathPoint > points( pathNodes.size() );
 
 
-//=======================================================================
-//function : makeEdgePathPoints
-//purpose  : auxiliary 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( FirstIsStart ) {
-    if ( aT1 > aT2 ) {
-      aPrms.reverse();
-    }
-  }
-  else {
-    if ( aT2 > aT1 ) {
-      aPrms.reverse();
-    }
-  }
-  // 3. Path Points
-  SMESH_MeshEditor_PathPoint aPP;
-  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( FirstIsStart ? aVec : -aVec );
-    aPP.SetPnt( aP3D );
-    aPP.SetTangent( aTgt );
-    aPP.SetParameter( aT );
-    LPP.push_back(aPP);
-  }
-  return EXTR_OK;
-}
+  if ( theAngleVariation )
+    linearAngleVariation( points.size()-1, theAngles );
+  if ( theScaleVariation )
+    linearScaleVariation( points.size()-1, theScales );
 
 
+  theAngles.push_front( 0 ); // for the 1st point that is not transformed
+  std::list<double>::iterator angle = theAngles.begin();
 
 
-//=======================================================================
-//function : makeExtrElements
-//purpose  : auxiliary for ExtrusionAlongTrack
-//=======================================================================
-SMESH_MeshEditor::Extrusion_Error
-SMESH_MeshEditor::makeExtrElements(TIDSortedElemSet                  theElemSets[2],
-                                   list<SMESH_MeshEditor_PathPoint>& fullList,
-                                   const bool                        theHasAngles,
-                                   list<double>&                     theAngles,
-                                   const bool                        theLinearVariation,
-                                   const bool                        theHasRefPoint,
-                                   const gp_Pnt&                     theRefPoint,
-                                   const bool                        theMakeGroups)
-{
-  const int aNbTP = fullList.size();
+  SMESHDS_Mesh* pathMeshDS = theTrackMesh->GetMeshDS();
 
 
-  // Angles
-  if( theHasAngles && !theAngles.empty() && theLinearVariation )
-    linearAngleVariation(aNbTP-1, theAngles);
+  std::map< int, double > edgeID2OriFactor; // orientation of EDGEs
+  std::map< int, double >::iterator id2factor;
+  SMESH_MesherHelper pathHelper( *theTrackMesh );
+  gp_Pnt p; gp_Vec tangent;
+  const double tol2 = gp::Resolution() * gp::Resolution();
 
 
-  // fill vector of path points with angles
-  vector<SMESH_MeshEditor_PathPoint> aPPs;
-  list<SMESH_MeshEditor_PathPoint>::iterator itPP = fullList.begin();
-  list<double>::iterator                 itAngles = theAngles.begin();
-  aPPs.push_back( *itPP++ );
-  for( ; itPP != fullList.end(); itPP++) {
-    aPPs.push_back( *itPP );
-    if ( theHasAngles && itAngles != theAngles.end() )
-      aPPs.back().SetAngle( *itAngles++ );
-  }
+  for ( size_t i = 0; i < pathNodes.size(); ++i )
+  {
+    ExtrusParam::PathPoint & point = points[ i ];
 
 
-  TNodeOfNodeListMap   mapNewNodes;
-  TElemOfVecOfNnlmiMap mapElemNewNodes;
-  TTElemOfElemListMap  newElemsMap;
-  TIDSortedElemSet::iterator itElem;
-  // source elements for each generated one
-  SMESH_SequenceOfElemPtr srcElems, srcNodes;
+    point.myPnt = SMESH_NodeXYZ( pathNodes[ i ]);
 
 
-  // 3. Center of rotation aV0
-  gp_Pnt aV0 = theRefPoint;
-  if ( !theHasRefPoint )
-  {
-    gp_XYZ aGC( 0.,0.,0. );
-    TIDSortedElemSet newNodes;
+    if ( angle != theAngles.end() )
+      point.myAngle = *angle++;
 
 
-    for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
+    tangent.SetCoord( 0,0,0 );
+    const int          shapeID = pathNodes[ i ]->GetShapeID();
+    const TopoDS_Shape&  shape = pathMeshDS->IndexToShape( shapeID );
+    TopAbs_ShapeEnum shapeType = shape.IsNull() ? TopAbs_SHAPE : shape.ShapeType();
+    switch ( shapeType )
     {
     {
-      TIDSortedElemSet& theElements = theElemSets[ is2ndSet ];
-      itElem = theElements.begin();
-      for ( ; itElem != theElements.end(); itElem++ )
+    case TopAbs_EDGE:
+    {
+      TopoDS_Edge edge = TopoDS::Edge( shape );
+      id2factor = edgeID2OriFactor.insert( std::make_pair( shapeID, 0 )).first;
+      if ( id2factor->second == 0 )
       {
       {
-        const SMDS_MeshElement* elem = *itElem;
-        SMDS_ElemIteratorPtr     itN = elem->nodesIterator();
-        while ( itN->more() ) {
-          const SMDS_MeshElement* node = itN->next();
-          if ( newNodes.insert( node ).second )
-            aGC += SMESH_NodeXYZ( node );
-        }
+        if ( i ) id2factor->second = getOriFactor( edge, pathNodes[i-1], pathNodes[i], pathHelper );
+        else     id2factor->second = getOriFactor( edge, pathNodes[i], pathNodes[i+1], pathHelper );
       }
       }
+      double u = pathHelper.GetNodeU( edge, pathNodes[i] ), u0, u1;
+      Handle(Geom_Curve) curve = BRep_Tool::Curve( edge, u0, u1 );
+      curve->D1( u, p, tangent );
+      tangent *= id2factor->second;
+      break;
     }
     }
-    aGC /= newNodes.size();
-    aV0.SetXYZ( aGC );
-  } // if (!theHasRefPoint) {
-
-  // 4. Processing the elements
-  SMESHDS_Mesh* aMesh = GetMeshDS();
-  list<const SMDS_MeshNode*> emptyList;
-
-  setElemsFirst( theElemSets );
-  for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
-  {
-    TIDSortedElemSet& theElements = theElemSets[ is2ndSet ];
-    for ( itElem = theElements.begin(); itElem != theElements.end(); itElem++ )
+    case TopAbs_VERTEX:
     {
     {
-      const SMDS_MeshElement* elem = *itElem;
-
-      vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
-      newNodesItVec.reserve( elem->NbNodes() );
-
-      // loop on elem nodes
-      int nodeIndex = -1;
-      SMDS_ElemIteratorPtr itN = elem->nodesIterator();
-      while ( itN->more() )
+      int nbEdges = 0;
+      PShapeIteratorPtr shapeIt = pathHelper.GetAncestors( shape, *theTrackMesh, TopAbs_EDGE );
+      while ( const TopoDS_Shape* edgePtr = shapeIt->next() )
       {
       {
-        ++nodeIndex;
-        // check if a node has been already processed
-        const SMDS_MeshNode* node = cast2Node( itN->next() );
-        TNodeOfNodeListMap::iterator nIt = mapNewNodes.insert( make_pair( node, emptyList )).first;
-        list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
-        if ( listNewNodes.empty() )
+        int edgeID = pathMeshDS->ShapeToIndex( *edgePtr );
+        for ( int di = -1; di <= 0; ++di )
         {
         {
-          // make new nodes
-          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 = SMESH_NodeXYZ( node );
-
-          const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0];
-          aP0x = aPP0.Pnt();
-          aDT0x= aPP0.Tangent();
-
-          for ( int 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 );
-
-            // translated point
-            aV1x = aV0x.Transformed( aTrsf );
-            aPN1 = aPN0.Transformed( aTrsf );
-
-            // rotation 1 [ T1,T0 ]
-            aAngleT1T0=-aDT1x.Angle( aDT0x );
-            if (fabs(aAngleT1T0) > aTolAng)
+          size_t j = i + di;
+          if ( j < pathEdges.size() && edgeID == pathEdges[ j ]->GetShapeID() )
+          {
+            TopoDS_Edge edge = TopoDS::Edge( *edgePtr );
+            id2factor = edgeID2OriFactor.insert( std::make_pair( edgeID, 0 )).first;
+            if ( id2factor->second == 0 )
             {
             {
-              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 );
+              if ( j < i )
+                id2factor->second = getOriFactor( edge, pathNodes[i-1], pathNodes[i], pathHelper );
+              else
+                id2factor->second = getOriFactor( edge, pathNodes[i], pathNodes[i+1], pathHelper );
             }
             }
-
-            // make new node
-            if ( elem->IsQuadratic() && !elem->IsMediumNode(node) )
+            double u = pathHelper.GetNodeU( edge, pathNodes[i] ), u0, u1;
+            Handle(Geom_Curve) curve = BRep_Tool::Curve( edge, u0, u1 );
+            gp_Vec du;
+            curve->D1( u, p, du );
+            double size2 = du.SquareMagnitude();
+            if ( du.SquareMagnitude() > tol2 )
             {
             {
-              // create additional node
-              gp_XYZ midP = 0.5 * ( aPN1.XYZ() + aPN0.XYZ() );
-              const SMDS_MeshNode* newNode = aMesh->AddNode( midP.X(), midP.Y(), midP.Z() );
-              myLastCreatedNodes.push_back(newNode);
-              srcNodes.push_back( node );
-              listNewNodes.push_back( newNode );
+              tangent += du.Divided( Sqrt( size2 )) * id2factor->second;
+              nbEdges++;
             }
             }
-            const SMDS_MeshNode* newNode = aMesh->AddNode( aPN1.X(), aPN1.Y(), aPN1.Z() );
-            myLastCreatedNodes.push_back(newNode);
-            srcNodes.push_back( node );
-            listNewNodes.push_back( newNode );
-
-            aPN0 = aPN1;
-            aP0x = aP1x;
-            aV0x = aV1x;
-            aDT0x = aDT1x;
+            break;
           }
         }
           }
         }
-        else if( elem->IsQuadratic() && !elem->IsMediumNode(node) )
+      }
+      if ( nbEdges > 0 )
+        break;
+    }
+    // fall through
+    default:
+    {
+      for ( int di = -1; di <= 1; di += 2 )
+      {
+        size_t j = i + di;
+        if ( j < pathNodes.size() )
         {
         {
-          // if current elem is quadratic and current node is not medium
-          // we have to check - may be it is needed to insert additional nodes
-          list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
-          if ((int) listNewNodes.size() == aNbTP-1 )
-          {
-            vector<const SMDS_MeshNode*> aNodes(2*(aNbTP-1));
-            gp_XYZ P(node->X(), node->Y(), node->Z());
-            list< const SMDS_MeshNode* >::iterator it = listNewNodes.begin();
-            int i;
-            for(i=0; i<aNbTP-1; i++) {
-              const SMDS_MeshNode* N = *it;
-              double x = ( N->X() + P.X() )/2.;
-              double y = ( N->Y() + P.Y() )/2.;
-              double z = ( N->Z() + P.Z() )/2.;
-              const SMDS_MeshNode* newN = aMesh->AddNode(x,y,z);
-              srcNodes.push_back( node );
-              myLastCreatedNodes.push_back(newN);
-              aNodes[2*i] = newN;
-              aNodes[2*i+1] = N;
-              P = gp_XYZ(N->X(),N->Y(),N->Z());
-            }
-            listNewNodes.clear();
-            for(i=0; i<2*(aNbTP-1); i++) {
-              listNewNodes.push_back(aNodes[i]);
-            }
-          }
+          gp_Vec dir( point.myPnt, SMESH_NodeXYZ( pathNodes[ j ]));
+          double size2 = dir.SquareMagnitude();
+          if ( size2 > tol2 )
+            tangent += dir.Divided( Sqrt( size2 )) * di;
         }
         }
-
-        newNodesItVec.push_back( nIt );
       }
       }
-
-      // make new elements
-      sweepElement( elem, newNodesItVec, newElemsMap[elem], aNbTP-1, srcElems );
     }
     }
-  }
+    } // switch ( shapeType )
 
 
-  makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElemSets[0], aNbTP-1, srcElems );
+    if ( tangent.SquareMagnitude() < tol2 )
+      return EXTR_CANT_GET_TANGENT;
 
 
-  if ( theMakeGroups )
-    generateGroups( srcNodes, srcElems, "extruded");
+    point.myTgt = tangent;
+
+  } // loop on pathNodes
+
+
+  ExtrusParam nodeMaker( points, theRefPoint, theScales, theMakeGroups );
+  TTElemOfElemListMap newElemsMap;
+
+  ExtrusionSweep( theElements, nodeMaker, newElemsMap );
 
   return EXTR_OK;
 }
 
 
   return EXTR_OK;
 }
 
-
 //=======================================================================
 //function : linearAngleVariation
 //purpose  : spread values over nbSteps
 //=======================================================================
 //function : linearAngleVariation
 //purpose  : spread values over nbSteps
@@ -6490,6 +6491,33 @@ void SMESH_MeshEditor::linearAngleVariation(const int     nbSteps,
   }
 }
 
   }
 }
 
+//=======================================================================
+//function : linearScaleVariation
+//purpose  : spread values over nbSteps 
+//=======================================================================
+
+void SMESH_MeshEditor::linearScaleVariation(const int          theNbSteps,
+                                            std::list<double>& theScales)
+{
+  int nbScales = theScales.size();
+  std::vector<double> myScales;
+  myScales.reserve( theNbSteps );
+  std::list<double>::const_iterator scale = theScales.begin();
+  double prevScale = 1.0;
+  for ( int iSc = 1; scale != theScales.end(); ++scale, ++iSc )
+  {
+    int      iStep = int( iSc / double( nbScales ) * theNbSteps + 0.5 );
+    int    stDelta = Max( 1, iStep - myScales.size());
+    double scDelta = ( *scale - prevScale ) / stDelta;
+    for ( int iStep = 0; iStep < stDelta; ++iStep )
+    {
+      myScales.push_back( prevScale + scDelta );
+      prevScale = myScales.back();
+    }
+    prevScale = *scale;
+  }
+  theScales.assign( myScales.begin(), myScales.end() );
+}
 
 //================================================================================
 /*!
 
 //================================================================================
 /*!
@@ -6732,7 +6760,7 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
  *  \param [in] theValue - offset value
  *  \param [out] theTgtMesh - a mesh to add offset elements to
  *  \param [in] theMakeGroups - to generate groups
  *  \param [in] theValue - offset value
  *  \param [out] theTgtMesh - a mesh to add offset elements to
  *  \param [in] theMakeGroups - to generate groups
- *  \return PGroupIDs - IDs of created groups
+ *  \return PGroupIDs - IDs of created groups. NULL means failure
  */
 //================================================================================
 
  */
 //================================================================================
 
@@ -6784,7 +6812,7 @@ SMESH_MeshEditor::PGroupIDs SMESH_MeshEditor::Offset( TIDSortedElemSet & theElem
 
   // copy offsetMesh to theTgtMesh
 
 
   // copy offsetMesh to theTgtMesh
 
-  int idShift = meshDS->MaxNodeID();
+  smIdType idShift = meshDS->MaxNodeID();
   for ( size_t i = 0; i < new2OldNodes.size(); ++i )
     if ( const SMDS_MeshNode* n = new2OldNodes[ i ].first )
     {
   for ( size_t i = 0; i < new2OldNodes.size(); ++i )
     if ( const SMDS_MeshNode* n = new2OldNodes[ i ].first )
     {
@@ -6819,6 +6847,8 @@ SMESH_MeshEditor::PGroupIDs SMESH_MeshEditor::Offset( TIDSortedElemSet & theElem
   PGroupIDs newGroupIDs;
   if ( theMakeGroups )
     newGroupIDs = generateGroups( srcNodes, srcElems, "offset", theTgtMesh, false );
   PGroupIDs newGroupIDs;
   if ( theMakeGroups )
     newGroupIDs = generateGroups( srcNodes, srcElems, "offset", theTgtMesh, false );
+  else
+    newGroupIDs.reset( new std::list< int > );
 
   return newGroupIDs;
 }
 
   return newGroupIDs;
 }
@@ -7152,7 +7182,7 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes,
 
   TNodeNodeMap nodeNodeMap; // node to replace - new node
   set<const SMDS_MeshElement*> elems; // all elements with changed nodes
 
   TNodeNodeMap nodeNodeMap; // node to replace - new node
   set<const SMDS_MeshElement*> elems; // all elements with changed nodes
-  list< int > rmElemIds, rmNodeIds;
+  list< smIdType > rmElemIds, rmNodeIds;
   vector< ElemFeatures > newElemDefs;
 
   // Fill nodeNodeMap and elems
   vector< ElemFeatures > newElemDefs;
 
   // Fill nodeNodeMap and elems
@@ -7258,9 +7288,19 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes,
 
     for ( size_t i = 0; i < newElemDefs.size(); ++i )
     {
 
     for ( size_t i = 0; i < newElemDefs.size(); ++i )
     {
-      if ( i > 0 || !mesh->ChangeElementNodes( elem,
-                                               & newElemDefs[i].myNodes[0],
-                                               newElemDefs[i].myNodes.size() ))
+      bool elemChanged = false;
+      if ( i == 0 )
+      {
+        if ( elem->GetGeomType() == SMDSGeom_POLYHEDRA )
+          elemChanged = mesh->ChangePolyhedronNodes( elem,
+                                                     newElemDefs[i].myNodes,
+                                                     newElemDefs[i].myPolyhedQuantities );
+        else
+          elemChanged = mesh->ChangeElementNodes( elem,
+                                                  & newElemDefs[i].myNodes[0],
+                                                  newElemDefs[i].myNodes.size() );
+      }
+      if ( i > 0 || !elemChanged )
       {
         if ( i == 0 )
         {
       {
         if ( i == 0 )
         {
@@ -7427,6 +7467,7 @@ bool SMESH_MeshEditor::applyMerge( const SMDS_MeshElement* elem,
         // each face has to be analyzed in order to check volume validity
         if ( const SMDS_MeshVolume* aPolyedre = SMDS_Mesh::DownCast< SMDS_MeshVolume >( elem ))
         {
         // each face has to be analyzed in order to check volume validity
         if ( const SMDS_MeshVolume* aPolyedre = SMDS_Mesh::DownCast< SMDS_MeshVolume >( elem ))
         {
+          toRemove = false;
           int nbFaces = aPolyedre->NbFaces();
 
           vector<const SMDS_MeshNode *>& poly_nodes = newElemDefs[0].myNodes;
           int nbFaces = aPolyedre->NbFaces();
 
           vector<const SMDS_MeshNode *>& poly_nodes = newElemDefs[0].myNodes;
@@ -7747,12 +7788,12 @@ bool SMESH_MeshEditor::applyMerge( const SMDS_MeshElement* elem,
 // purpose : allow comparing elements basing on their nodes
 // ========================================================
 
 // purpose : allow comparing elements basing on their nodes
 // ========================================================
 
-class ComparableElement : public boost::container::flat_set< int >
+class ComparableElement : public boost::container::flat_set< smIdType >
 {
 {
-  typedef boost::container::flat_set< int >  int_set;
+  typedef boost::container::flat_set< smIdType >  int_set;
 
   const SMDS_MeshElement* myElem;
 
   const SMDS_MeshElement* myElem;
-  int                     mySumID;
+  smIdType                mySumID;
   mutable int             myGroupID;
 
 public:
   mutable int             myGroupID;
 
 public:
@@ -7763,7 +7804,7 @@ public:
     this->reserve( theElem->NbNodes() );
     for ( SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator(); nodeIt->more(); )
     {
     this->reserve( theElem->NbNodes() );
     for ( SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator(); nodeIt->more(); )
     {
-      int id = nodeIt->next()->GetID();
+      smIdType id = nodeIt->next()->GetID();
       mySumID += id;
       this->insert( id );
     }
       mySumID += id;
       this->insert( id );
     }
@@ -7775,9 +7816,10 @@ public:
   //int& GroupID() const { return const_cast< int& >( myGroupID ); }
 
   ComparableElement( const ComparableElement& theSource ) // move copy
   //int& GroupID() const { return const_cast< int& >( myGroupID ); }
 
   ComparableElement( const ComparableElement& theSource ) // move copy
+    : int_set()
   {
     ComparableElement& src = const_cast< ComparableElement& >( theSource );
   {
     ComparableElement& src = const_cast< ComparableElement& >( theSource );
-    (int_set&) (*this ) = boost::move( src );
+    (int_set&) (*this ) = std::move( src );
     myElem    = src.myElem;
     mySumID   = src.mySumID;
     myGroupID = src.myGroupID;
     myElem    = src.myElem;
     mySumID   = src.mySumID;
     myGroupID = src.myGroupID;
@@ -7785,7 +7827,7 @@ public:
 
   static int HashCode(const ComparableElement& se, int limit )
   {
 
   static int HashCode(const ComparableElement& se, int limit )
   {
-    return ::HashCode( se.mySumID, limit );
+    return ::HashCode( FromSmIdType<int>(se.mySumID), limit );
   }
   static Standard_Boolean IsEqual(const ComparableElement& se1, const ComparableElement& se2 )
   {
   }
   static Standard_Boolean IsEqual(const ComparableElement& se1, const ComparableElement& se2 )
   {
@@ -7810,7 +7852,7 @@ void SMESH_MeshEditor::FindEqualElements( TIDSortedElemSet &        theElements,
   else                       elemIt = SMESHUtils::elemSetIterator( theElements );
 
   typedef NCollection_Map< ComparableElement, ComparableElement > TMapOfElements;
   else                       elemIt = SMESHUtils::elemSetIterator( theElements );
 
   typedef NCollection_Map< ComparableElement, ComparableElement > TMapOfElements;
-  typedef std::list<int>                                          TGroupOfElems;
+  typedef std::list<smIdType>                                     TGroupOfElems;
   TMapOfElements               mapOfElements;
   std::vector< TGroupOfElems > arrayOfGroups;
   TGroupOfElems                groupOfElems;
   TMapOfElements               mapOfElements;
   std::vector< TGroupOfElems > arrayOfGroups;
   TGroupOfElems                groupOfElems;
@@ -7818,6 +7860,8 @@ void SMESH_MeshEditor::FindEqualElements( TIDSortedElemSet &        theElements,
   while ( elemIt->more() )
   {
     const SMDS_MeshElement* curElem = elemIt->next();
   while ( elemIt->more() )
   {
     const SMDS_MeshElement* curElem = elemIt->next();
+    if ( curElem->IsNull() )
+      continue;
     ComparableElement      compElem = curElem;
     // check uniqueness
     const ComparableElement& elemInSet = mapOfElements.Added( compElem );
     ComparableElement      compElem = curElem;
     // check uniqueness
     const ComparableElement& elemInSet = mapOfElements.Added( compElem );
@@ -7854,7 +7898,7 @@ void SMESH_MeshEditor::MergeElements(TListOfListOfElementsID & theGroupsOfElemen
 {
   ClearLastCreated();
 
 {
   ClearLastCreated();
 
-  typedef list<int> TListOfIDs;
+  typedef list<smIdType> TListOfIDs;
   TListOfIDs rmElemIds; // IDs of elems to remove
 
   SMESHDS_Mesh* aMesh = GetMeshDS();
   TListOfIDs rmElemIds; // IDs of elems to remove
 
   SMESHDS_Mesh* aMesh = GetMeshDS();
@@ -7953,35 +7997,37 @@ bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode*             theFirst
   theNodes.push_back( theSecondNode );
 
   const SMDS_MeshNode *nIgnore = theFirstNode, *nStart = theSecondNode;
   theNodes.push_back( theSecondNode );
 
   const SMDS_MeshNode *nIgnore = theFirstNode, *nStart = theSecondNode;
-  TIDSortedElemSet foundElems;
+  //TIDSortedElemSet foundElems;
   bool needTheLast = ( theLastNode != 0 );
 
   bool needTheLast = ( theLastNode != 0 );
 
+  vector<const SMDS_MeshNode*> nodes;
+  
   while ( nStart != theLastNode ) {
     if ( nStart == theFirstNode )
       return !needTheLast;
 
   while ( nStart != theLastNode ) {
     if ( nStart == theFirstNode )
       return !needTheLast;
 
-    // find all free border faces sharing form nStart
+    // find all free border faces sharing nStart
 
     list< const SMDS_MeshElement* > curElemList;
     list< const SMDS_MeshNode* >    nStartList;
     SMDS_ElemIteratorPtr invElemIt = nStart->GetInverseElementIterator(SMDSAbs_Face);
     while ( invElemIt->more() ) {
       const SMDS_MeshElement* e = invElemIt->next();
 
     list< const SMDS_MeshElement* > curElemList;
     list< const SMDS_MeshNode* >    nStartList;
     SMDS_ElemIteratorPtr invElemIt = nStart->GetInverseElementIterator(SMDSAbs_Face);
     while ( invElemIt->more() ) {
       const SMDS_MeshElement* e = invElemIt->next();
-      if ( e == curElem || foundElems.insert( e ).second ) {
+      //if ( e == curElem || foundElems.insert( e ).second ) // e can encounter twice in border
+      {
         // get nodes
         // get nodes
-        int iNode = 0, nbNodes = e->NbNodes();
-        vector<const SMDS_MeshNode*> nodes( nbNodes+1 );
         nodes.assign( SMDS_MeshElement::iterator( e->interlacedNodesIterator() ),
                       SMDS_MeshElement::iterator() );
         nodes.push_back( nodes[ 0 ]);
 
         // check 2 links
         nodes.assign( SMDS_MeshElement::iterator( e->interlacedNodesIterator() ),
                       SMDS_MeshElement::iterator() );
         nodes.push_back( nodes[ 0 ]);
 
         // check 2 links
+        int iNode = 0, nbNodes = nodes.size() - 1;
         for ( iNode = 0; iNode < nbNodes; iNode++ )
         for ( iNode = 0; iNode < nbNodes; iNode++ )
-          if (((nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) ||
-               (nodes[ iNode + 1] == nStart && nodes[ iNode ] != nIgnore )) &&
-              ControlFreeBorder( &nodes[ iNode ], e->GetID() ))
+          if ((( nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) ||
+               ( nodes[ iNode + 1] == nStart && nodes[ iNode ] != nIgnore )) &&
+              ( ControlFreeBorder( &nodes[ iNode ], e->GetID() )))
           {
           {
-            nStartList.push_back( nodes[ iNode + ( nodes[ iNode ] == nStart ? 1 : 0 )]);
+            nStartList.push_back( nodes[ iNode + ( nodes[ iNode ] == nStart )]);
             curElemList.push_back( e );
           }
       }
             curElemList.push_back( e );
           }
       }
@@ -8033,7 +8079,7 @@ bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode*             theFirst
         else if ( !contNodes[0].empty() && !contNodes[1].empty() ) {
           // choice: clear a worse one
           int iLongest = ( contNodes[0].size() < contNodes[1].size() ? 1 : 0 );
         else if ( !contNodes[0].empty() && !contNodes[1].empty() ) {
           // choice: clear a worse one
           int iLongest = ( contNodes[0].size() < contNodes[1].size() ? 1 : 0 );
-          int iWorse = ( needTheLast ? 1 - iLongest : iLongest );
+          int   iWorse = ( needTheLast ? 1 - iLongest : iLongest );
           contNodes[ iWorse ].clear();
           contFaces[ iWorse ].clear();
         }
           contNodes[ iWorse ].clear();
           contFaces[ iWorse ].clear();
         }
@@ -8044,13 +8090,11 @@ bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode*             theFirst
       // push_back the best free border
       cNL = & contNodes[ contNodes[0].empty() ? 1 : 0 ];
       cFL = & contFaces[ contFaces[0].empty() ? 1 : 0 ];
       // push_back the best free border
       cNL = & contNodes[ contNodes[0].empty() ? 1 : 0 ];
       cFL = & contFaces[ contFaces[0].empty() ? 1 : 0 ];
-      theNodes.pop_back(); // remove nIgnore
+      //theNodes.pop_back(); // remove nIgnore
       theNodes.pop_back(); // remove nStart
       theNodes.pop_back(); // remove nStart
-      theFaces.pop_back(); // remove curElem
-      list< const SMDS_MeshNode* >::iterator nIt = cNL->begin();
-      list< const SMDS_MeshElement* >::iterator fIt = cFL->begin();
-      for ( ; nIt != cNL->end(); nIt++ ) theNodes.push_back( *nIt );
-      for ( ; fIt != cFL->end(); fIt++ ) theFaces.push_back( *fIt );
+      //theFaces.pop_back(); // remove curElem
+      theNodes.splice( theNodes.end(), *cNL );
+      theFaces.splice( theFaces.end(), *cFL );
       return true;
 
     } // several continuations found
       return true;
 
     } // several continuations found
@@ -8391,6 +8435,10 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
     nIt[0] = nSide[0].begin(); eIt[0] = eSide[0].begin();
     nIt[1] = nSide[1].begin(); eIt[1] = eSide[1].begin();
 
     nIt[0] = nSide[0].begin(); eIt[0] = eSide[0].begin();
     nIt[1] = nSide[1].begin(); eIt[1] = eSide[1].begin();
 
+    // element can be split while iterating on border if it has two edges in the border
+    std::map< const SMDS_MeshElement* , const SMDS_MeshElement* > elemReplaceMap;
+    std::map< const SMDS_MeshElement* , const SMDS_MeshElement* >::iterator elemReplaceMapIt;
+
     TElemOfNodeListMap insertMap;
     TElemOfNodeListMap::iterator insertMapIt;
     // insertMap is
     TElemOfNodeListMap insertMap;
     TElemOfNodeListMap::iterator insertMapIt;
     // insertMap is
@@ -8438,12 +8486,15 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
         const SMDS_MeshNode*    nIns = *nIt [ 1 - intoBord ];
         if ( intoBord == 1 ) {
           // move node of the border to be on a link of elem of the side
         const SMDS_MeshNode*    nIns = *nIt [ 1 - intoBord ];
         if ( intoBord == 1 ) {
           // move node of the border to be on a link of elem of the side
-          gp_XYZ p1 (n1->X(), n1->Y(), n1->Z());
-          gp_XYZ p2 (n2->X(), n2->Y(), n2->Z());
+          SMESH_NodeXYZ p1( n1 ), p2( n2 );
           double ratio = du / ( param[ 1 ][ i[1] ] - param[ 1 ][ i[1]-1 ]);
           gp_XYZ p = p2 * ( 1 - ratio ) + p1 * ratio;
           GetMeshDS()->MoveNode( nIns, p.X(), p.Y(), p.Z() );
         }
           double ratio = du / ( param[ 1 ][ i[1] ] - param[ 1 ][ i[1]-1 ]);
           gp_XYZ p = p2 * ( 1 - ratio ) + p1 * ratio;
           GetMeshDS()->MoveNode( nIns, p.X(), p.Y(), p.Z() );
         }
+        elemReplaceMapIt = elemReplaceMap.find( elem );
+        if ( elemReplaceMapIt != elemReplaceMap.end() )
+          elem = elemReplaceMapIt->second;
+
         insertMapIt = insertMap.find( elem );
         bool  notFound = ( insertMapIt == insertMap.end() );
         bool otherLink = ( !notFound && (*insertMapIt).second.front() != n1 );
         insertMapIt = insertMap.find( elem );
         bool  notFound = ( insertMapIt == insertMap.end() );
         bool otherLink = ( !notFound && (*insertMapIt).second.front() != n1 );
@@ -8466,8 +8517,10 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
             UpdateVolumes(n12, n22, nodeList);
           }
           // 3. find an element appeared on n1 and n2 after the insertion
             UpdateVolumes(n12, n22, nodeList);
           }
           // 3. find an element appeared on n1 and n2 after the insertion
-          insertMap.erase( elem );
-          elem = findAdjacentFace( n1, n2, 0 );
+          insertMap.erase( insertMapIt );
+          const SMDS_MeshElement* elem2 = findAdjacentFace( n1, n2, 0 );
+          elemReplaceMap.insert( std::make_pair( elem, elem2 ));
+          elem = elem2;
         }
         if ( notFound || otherLink ) {
           // add element and nodes of the side into the insertMap
         }
         if ( notFound || otherLink ) {
           // add element and nodes of the side into the insertMap
@@ -8560,8 +8613,8 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
     TListOfListOfElementsID::iterator itGroups = equalGroups.begin();
     for ( ; itGroups != equalGroups.end(); ++itGroups )
     {
     TListOfListOfElementsID::iterator itGroups = equalGroups.begin();
     for ( ; itGroups != equalGroups.end(); ++itGroups )
     {
-      list< int >& group = *itGroups;
-      list< int >::iterator id = group.begin();
+      list< smIdType >& group = *itGroups;
+      list< smIdType >::iterator id = group.begin();
       for ( ++id; id != group.end(); ++id )
         if ( const SMDS_MeshElement* seg = GetMeshDS()->FindElement( *id ))
           segments.erase( seg );
       for ( ++id; id != group.end(); ++id )
         if ( const SMDS_MeshElement* seg = GetMeshDS()->FindElement( *id ))
           segments.erase( seg );
@@ -8993,12 +9046,12 @@ namespace
  */
 //=======================================================================
 
  */
 //=======================================================================
 
-int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
-                                             SMESH_MesherHelper& theHelper,
-                                             const bool          theForce3d)
+smIdType SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
+                                                  SMESH_MesherHelper& theHelper,
+                                                  const bool          theForce3d)
 {
   //MESSAGE("convertElemToQuadratic");
 {
   //MESSAGE("convertElemToQuadratic");
-  int nbElem = 0;
+  smIdType nbElem = 0;
   if( !theSm ) return nbElem;
 
   vector<int> nbNodeInFaces;
   if( !theSm ) return nbElem;
 
   vector<int> nbNodeInFaces;
@@ -9051,7 +9104,7 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
     }
     // get elem data needed to re-create it
     //
     }
     // get elem data needed to re-create it
     //
-    const int id      = elem->GetID();
+    const smIdType id = elem->GetID();
     const int nbNodes = elem->NbCornerNodes();
     nodes.assign(elem->begin_nodes(), elem->end_nodes());
     if ( aGeomType == SMDSEntity_Polyhedra )
     const int nbNodes = elem->NbCornerNodes();
     nodes.assign(elem->begin_nodes(), elem->end_nodes());
     if ( aGeomType == SMDSEntity_Polyhedra )
@@ -9146,7 +9199,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT
   aHelper.ToFixNodeParameters( true );
 
   // convert elements assigned to sub-meshes
   aHelper.ToFixNodeParameters( true );
 
   // convert elements assigned to sub-meshes
-  int nbCheckedElems = 0;
+  smIdType nbCheckedElems = 0;
   if ( myMesh->HasShapeToMesh() )
   {
     if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
   if ( myMesh->HasShapeToMesh() )
   {
     if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
@@ -9163,7 +9216,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT
   }
 
   // convert elements NOT assigned to sub-meshes
   }
 
   // convert elements NOT assigned to sub-meshes
-  int totalNbElems = meshDS->NbEdges() + meshDS->NbFaces() + meshDS->NbVolumes();
+  smIdType totalNbElems = meshDS->NbEdges() + meshDS->NbFaces() + meshDS->NbVolumes();
   if ( nbCheckedElems < totalNbElems ) // not all elements are in sub-meshes
   {
     aHelper.SetElementsOnShape(false);
   if ( nbCheckedElems < totalNbElems ) // not all elements are in sub-meshes
   {
     aHelper.SetElementsOnShape(false);
@@ -9176,7 +9229,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT
       const SMDS_MeshEdge* edge = aEdgeItr->next();
       if ( !edge->IsQuadratic() )
       {
       const SMDS_MeshEdge* edge = aEdgeItr->next();
       if ( !edge->IsQuadratic() )
       {
-        int                  id = edge->GetID();
+        smIdType                  id = edge->GetID();
         const SMDS_MeshNode* n1 = edge->GetNode(0);
         const SMDS_MeshNode* n2 = edge->GetNode(1);
 
         const SMDS_MeshNode* n1 = edge->GetNode(0);
         const SMDS_MeshNode* n2 = edge->GetNode(1);
 
@@ -9217,7 +9270,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT
       if ( alreadyOK )
         continue;
 
       if ( alreadyOK )
         continue;
 
-      const int id = face->GetID();
+      const smIdType id = face->GetID();
       vector<const SMDS_MeshNode *> nodes ( face->begin_nodes(), face->end_nodes());
 
       meshDS->RemoveFreeElement(face, smDS, /*fromGroups=*/false);
       vector<const SMDS_MeshNode *> nodes ( face->begin_nodes(), face->end_nodes());
 
       meshDS->RemoveFreeElement(face, smDS, /*fromGroups=*/false);
@@ -9273,7 +9326,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT
           continue;
         }
       }
           continue;
         }
       }
-      const int id = volume->GetID();
+      const smIdType id = volume->GetID();
       vector<const SMDS_MeshNode *> nodes (volume->begin_nodes(), volume->end_nodes());
       if ( type == SMDSEntity_Polyhedra )
         nbNodeInFaces = static_cast<const SMDS_MeshVolume* >(volume)->GetQuantities();
       vector<const SMDS_MeshNode *> nodes (volume->begin_nodes(), volume->end_nodes());
       if ( type == SMDSEntity_Polyhedra )
         nbNodeInFaces = static_cast<const SMDS_MeshVolume* >(volume)->GetQuantities();
@@ -9453,7 +9506,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool        theForce3d,
     if ( alreadyOK ) continue;
 
     const SMDSAbs_ElementType type = elem->GetType();
     if ( alreadyOK ) continue;
 
     const SMDSAbs_ElementType type = elem->GetType();
-    const int                   id = elem->GetID();
+    const smIdType              id = elem->GetID();
     const int              nbNodes = elem->NbCornerNodes();
     vector<const SMDS_MeshNode *> nodes ( elem->begin_nodes(), elem->end_nodes());
 
     const int              nbNodes = elem->NbCornerNodes();
     vector<const SMDS_MeshNode *> nodes ( elem->begin_nodes(), elem->end_nodes());
 
@@ -9514,15 +9567,15 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool        theForce3d,
 //=======================================================================
 /*!
  * \brief Convert quadratic elements to linear ones and remove quadratic nodes
 //=======================================================================
 /*!
  * \brief Convert quadratic elements to linear ones and remove quadratic nodes
- * \return int - nb of checked elements
+ * \return smIdType - nb of checked elements
  */
 //=======================================================================
 
  */
 //=======================================================================
 
-int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh *    theSm,
-                                     SMDS_ElemIteratorPtr theItr,
-                                     const int            theShapeID)
+smIdType SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh *    theSm,
+                                          SMDS_ElemIteratorPtr theItr,
+                                          const int            /*theShapeID*/)
 {
 {
-  int nbElem = 0;
+  smIdType nbElem = 0;
   SMESHDS_Mesh* meshDS = GetMeshDS();
   ElemFeatures elemType;
   vector<const SMDS_MeshNode *> nodes;
   SMESHDS_Mesh* meshDS = GetMeshDS();
   ElemFeatures elemType;
   vector<const SMDS_MeshNode *> nodes;
@@ -9567,7 +9620,7 @@ int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh *    theSm,
 
 bool SMESH_MeshEditor::ConvertFromQuadratic()
 {
 
 bool SMESH_MeshEditor::ConvertFromQuadratic()
 {
-  int nbCheckedElems = 0;
+  smIdType nbCheckedElems = 0;
   if ( myMesh->HasShapeToMesh() )
   {
     if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
   if ( myMesh->HasShapeToMesh() )
   {
     if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
@@ -9581,7 +9634,7 @@ bool SMESH_MeshEditor::ConvertFromQuadratic()
     }
   }
 
     }
   }
 
-  int totalNbElems =
+  smIdType totalNbElems =
     GetMeshDS()->NbEdges() + GetMeshDS()->NbFaces() + GetMeshDS()->NbVolumes();
   if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes
   {
     GetMeshDS()->NbEdges() + GetMeshDS()->NbFaces() + GetMeshDS()->NbVolumes();
   if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes
   {
@@ -9620,7 +9673,7 @@ void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements)
   if ( theElements.empty() ) return;
 
   // collect IDs of medium nodes of theElements; some of these nodes will be removed
   if ( theElements.empty() ) return;
 
   // collect IDs of medium nodes of theElements; some of these nodes will be removed
-  set<int> mediumNodeIDs;
+  set<smIdType> mediumNodeIDs;
   TIDSortedElemSet::iterator eIt = theElements.begin();
   for ( ; eIt != theElements.end(); ++eIt )
   {
   TIDSortedElemSet::iterator eIt = theElements.begin();
   for ( ; eIt != theElements.end(); ++eIt )
   {
@@ -9639,7 +9692,7 @@ void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements)
 
   // get remaining medium nodes
   TIDSortedNodeSet mediumNodes;
 
   // get remaining medium nodes
   TIDSortedNodeSet mediumNodes;
-  set<int>::iterator nIdsIt = mediumNodeIDs.begin();
+  set<smIdType>::iterator nIdsIt = mediumNodeIDs.begin();
   for ( ; nIdsIt != mediumNodeIDs.end(); ++nIdsIt )
     if ( const SMDS_MeshNode* n = GetMeshDS()->FindNode( *nIdsIt ))
       mediumNodes.insert( mediumNodes.end(), n );
   for ( ; nIdsIt != mediumNodeIDs.end(); ++nIdsIt )
     if ( const SMDS_MeshNode* n = GetMeshDS()->FindNode( *nIdsIt ))
       mediumNodes.insert( mediumNodes.end(), n );
@@ -9667,7 +9720,7 @@ void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements)
             const SMDS_MeshElement* eComplex = invIt2->next();
             if ( eComplex->IsQuadratic() && !allMediumNodesIn( eComplex, mediumNodes))
             {
             const SMDS_MeshElement* eComplex = invIt2->next();
             if ( eComplex->IsQuadratic() && !allMediumNodesIn( eComplex, mediumNodes))
             {
-              int nbCommonNodes = SMESH_MeshAlgos::GetCommonNodes( e, eComplex ).size();
+              int nbCommonNodes = SMESH_MeshAlgos::NbCommonNodes( e, eComplex );
               if ( nbCommonNodes == e->NbNodes())
               {
                 complexFound = true;
               if ( nbCommonNodes == e->NbNodes())
               {
                 complexFound = true;
@@ -10128,7 +10181,7 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
   if ( aResult != SEW_OK)
     return aResult;
 
   if ( aResult != SEW_OK)
     return aResult;
 
-  list< int > nodeIDsToRemove;
+  list< smIdType > nodeIDsToRemove;
   vector< const SMDS_MeshNode*> nodes;
   ElemFeatures elemType;
 
   vector< const SMDS_MeshNode*> nodes;
   ElemFeatures elemType;
 
@@ -10209,10 +10262,10 @@ SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
   set<const SMDS_MeshElement*> * faceSetPtr[] = { &theSide1, &theSide2 };
 
   nReplaceMap.clear();
   set<const SMDS_MeshElement*> * faceSetPtr[] = { &theSide1, &theSide2 };
 
   nReplaceMap.clear();
-  if ( theFirstNode1 != theFirstNode2 )
-    nReplaceMap.insert( make_pair( theFirstNode1, theFirstNode2 ));
-  if ( theSecondNode1 != theSecondNode2 )
-    nReplaceMap.insert( make_pair( theSecondNode1, theSecondNode2 ));
+  //if ( theFirstNode1 != theFirstNode2 )
+  nReplaceMap.insert( make_pair( theFirstNode1, theFirstNode2 ));
+  //if ( theSecondNode1 != theSecondNode2 )
+  nReplaceMap.insert( make_pair( theSecondNode1, theSecondNode2 ));
 
   set< SMESH_TLink > linkSet; // set of nodes where order of nodes is ignored
   linkSet.insert( SMESH_TLink( theFirstNode1, theSecondNode1 ));
 
   set< SMESH_TLink > linkSet; // set of nodes where order of nodes is ignored
   linkSet.insert( SMESH_TLink( theFirstNode1, theSecondNode1 ));
@@ -10502,7 +10555,7 @@ namespace // automatically find theAffectedElems for DoubleNodes()
         if ( SMESH_MeshAlgos::FaceNormal( _elems[1], norm ))
           avgNorm += norm;
 
         if ( SMESH_MeshAlgos::FaceNormal( _elems[1], norm ))
           avgNorm += norm;
 
-        gp_XYZ bordDir( SMESH_NodeXYZ( _nodes[0] ) - SMESH_NodeXYZ( _nodes[1] ));
+        gp_XYZ bordDir( SMESH_NodeXYZ( this->_nodes[0] ) - SMESH_NodeXYZ( this->_nodes[1] ));
         norm = bordDir ^ avgNorm;
       }
       else
         norm = bordDir ^ avgNorm;
       }
       else
@@ -10528,7 +10581,8 @@ namespace // automatically find theAffectedElems for DoubleNodes()
       if ( maxX < 0 )
       {
         _elems[0]->setIsMarked( false );
       if ( maxX < 0 )
       {
         _elems[0]->setIsMarked( false );
-        _elems[1]->setIsMarked( true );
+        if ( _elems[1] )
+          _elems[1]->setIsMarked( true );
       }
     }
 
       }
     }
 
@@ -10700,7 +10754,11 @@ namespace // automatically find theAffectedElems for DoubleNodes()
     {
       fissure.reserve( theElemsOrNodes.size() );
       for ( ; elIt != theElemsOrNodes.end(); ++elIt )
     {
       fissure.reserve( theElemsOrNodes.size() );
       for ( ; elIt != theElemsOrNodes.end(); ++elIt )
+      {
         fissure.push_back( std::move( FissureBorder( *elIt, elemsByFacet )));
         fissure.push_back( std::move( FissureBorder( *elIt, elemsByFacet )));
+        if ( !fissure.back()._elems[1] )
+          fissure.pop_back();
+      }
     }
     if ( fissure.empty() )
       return;
     }
     if ( fissure.empty() )
       return;
@@ -11063,7 +11121,23 @@ bool SMESH_MeshEditor::doubleNodes(SMESHDS_Mesh*           theMeshDS,
     if ( theIsDoubleElem )
       AddElement( newNodes, elemType.Init( anElem, /*basicOnly=*/false ));
     else
     if ( theIsDoubleElem )
       AddElement( newNodes, elemType.Init( anElem, /*basicOnly=*/false ));
     else
-      theMeshDS->ChangeElementNodes( anElem, &newNodes[ 0 ], newNodes.size() );
+    {
+      // change element nodes
+      const SMDSAbs_EntityType geomType = anElem->GetEntityType();
+      if ( geomType == SMDSEntity_Polyhedra )
+      {
+        // special treatment for polyhedron
+        const SMDS_MeshVolume* aPolyhedron = SMDS_Mesh::DownCast< SMDS_MeshVolume >( anElem );
+        if (!aPolyhedron) {
+          MESSAGE("Warning: bad volumic element");
+          return false;
+        }
+        theMeshDS->ChangePolyhedronNodes( anElem, newNodes, aPolyhedron->GetQuantities() );
+      }
+      else
+        // standard entity type
+        theMeshDS->ChangeElementNodes( anElem, &newNodes[ 0 ], newNodes.size() );
+    }
 
     res = true;
   }
 
     res = true;
   }
@@ -11207,7 +11281,7 @@ namespace {
   \brief Identify the elements that will be affected by node duplication (actual duplication is not performed).
   This method is the first step of DoubleNodeElemGroupsInRegion.
   \param theElems - list of groups of elements (edges or faces) to be replicated
   \brief Identify the elements that will be affected by node duplication (actual duplication is not performed).
   This method is the first step of DoubleNodeElemGroupsInRegion.
   \param theElems - list of groups of elements (edges or faces) to be replicated
-  \param theNodesNot - list of groups of nodes not to replicated
+  \param theNodesNot - list of groups of nodes not to replicate
   \param theShape - shape to detect affected elements (element which geometric center
          located on or inside shape). If the shape is null, detection is done on faces orientations
          (select elements with a gravity center on the side given by faces normals).
   \param theShape - shape to detect affected elements (element which geometric center
          located on or inside shape). If the shape is null, detection is done on faces orientations
          (select elements with a gravity center on the side given by faces normals).
@@ -11349,7 +11423,7 @@ double SMESH_MeshEditor::OrientedAngle(const gp_Pnt& p0, const gp_Pnt& p1, const
   try {
     return n2.AngleWithRef(n1, vref);
   }
   try {
     return n2.AngleWithRef(n1, vref);
   }
-  catch ( Standard_Failure ) {
+  catch ( Standard_Failure& ) {
   }
   return Max( v1.Magnitude(), v2.Magnitude() );
 }
   }
   return Max( v1.Magnitude(), v2.Magnitude() );
 }
@@ -11389,16 +11463,9 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   //     build the list of nodes shared by 2 or more domains, with their domain indexes
 
   std::map<DownIdType, std::map<int,int>, DownIdCompare> faceDomains; // face --> (id domain --> id volume)
   //     build the list of nodes shared by 2 or more domains, with their domain indexes
 
   std::map<DownIdType, std::map<int,int>, DownIdCompare> faceDomains; // face --> (id domain --> id volume)
-  std::map<int,int>celldom; // cell vtkId --> domain
+  std::map<int,int> celldom; // cell vtkId --> domain
   std::map<DownIdType, std::map<int,int>, DownIdCompare> cellDomains;  // oldNode --> (id domain --> id cell)
   std::map<int, std::map<int,int> > nodeDomains; // oldId -->  (domainId --> newId)
   std::map<DownIdType, std::map<int,int>, DownIdCompare> cellDomains;  // oldNode --> (id domain --> id cell)
   std::map<int, std::map<int,int> > nodeDomains; // oldId -->  (domainId --> newId)
-  faceDomains.clear();
-  celldom.clear();
-  cellDomains.clear();
-  nodeDomains.clear();
-  std::map<int,int> emptyMap;
-  std::set<int> emptySet;
-  emptyMap.clear();
 
   //MESSAGE(".. Number of domains :"<<theElems.size());
 
 
   //MESSAGE(".. Number of domains :"<<theElems.size());
 
@@ -11448,7 +11515,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       const SMDS_MeshElement* anElem = *elemItr;
       if (!anElem)
         continue;
       const SMDS_MeshElement* anElem = *elemItr;
       if (!anElem)
         continue;
-      int vtkId = anElem->GetVtkID();
+      vtkIdType vtkId = anElem->GetVtkID();
       //MESSAGE("  vtkId " << vtkId << " smdsId " << anElem->GetID());
       int neighborsVtkIds[NBMAXNEIGHBORS];
       int downIds[NBMAXNEIGHBORS];
       //MESSAGE("  vtkId " << vtkId << " smdsId " << anElem->GetID());
       int neighborsVtkIds[NBMAXNEIGHBORS];
       int downIds[NBMAXNEIGHBORS];
@@ -11456,7 +11523,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId);
       for (int n = 0; n < nbNeighbors; n++)
       {
       int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId);
       for (int n = 0; n < nbNeighbors; n++)
       {
-        int smdsId = meshDS->FromVtkToSmds(neighborsVtkIds[n]);
+        smIdType smdsId = meshDS->FromVtkToSmds(neighborsVtkIds[n]);
         const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
         if (elem && ! domain.count(elem)) // neighbor is in another domain : face is shared
         {
         const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
         if (elem && ! domain.count(elem)) // neighbor is in another domain : face is shared
         {
@@ -11508,14 +11575,13 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       DownIdType face = itface->first;
       //MESSAGE(" --- face " << face.cellId);
       std::set<int> oldNodes;
       DownIdType face = itface->first;
       //MESSAGE(" --- face " << face.cellId);
       std::set<int> oldNodes;
-      oldNodes.clear();
       grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
       std::set<int>::iterator itn = oldNodes.begin();
       for (; itn != oldNodes.end(); ++itn)
       {
         int oldId = *itn;
         //MESSAGE("     node " << oldId);
       grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
       std::set<int>::iterator itn = oldNodes.begin();
       for (; itn != oldNodes.end(); ++itn)
       {
         int oldId = *itn;
         //MESSAGE("     node " << oldId);
-        vtkCellLinks::Link l = grid->GetCellLinks()->GetLink(oldId);
+        vtkCellLinks::Link l = (static_cast <vtkCellLinks *>(grid->GetCellLinks()))->GetLink(oldId);
         for (int i=0; i<l.ncells; i++)
         {
           int vtkId = l.cells[i];
         for (int i=0; i<l.ncells; i++)
         {
           int vtkId = l.cells[i];
@@ -11563,7 +11629,6 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       DownIdType face = itface->first;
       //MESSAGE(" --- face " << face.cellId);
       std::set<int> oldNodes;
       DownIdType face = itface->first;
       //MESSAGE(" --- face " << face.cellId);
       std::set<int> oldNodes;
-      oldNodes.clear();
       grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
       std::set<int>::iterator itn = oldNodes.begin();
       for (; itn != oldNodes.end(); ++itn)
       grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
       std::set<int>::iterator itn = oldNodes.begin();
       for (; itn != oldNodes.end(); ++itn)
@@ -11624,7 +11689,6 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       DownIdType face = itface->first;
       //MESSAGE(" --- face " << face.cellId);
       std::set<int> oldNodes;
       DownIdType face = itface->first;
       //MESSAGE(" --- face " << face.cellId);
       std::set<int> oldNodes;
-      oldNodes.clear();
       grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
       int nbMultipleNodes = 0;
       std::set<int>::iterator itn = oldNodes.begin();
       grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
       int nbMultipleNodes = 0;
       std::set<int>::iterator itn = oldNodes.begin();
@@ -11685,7 +11749,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
                   const TIDSortedElemSet& domain = (idom == iRestDom) ? theRestDomElems : theElems[idom];
                   for ( int ivol = 0; ivol < nbvol; ivol++ )
                   {
                   const TIDSortedElemSet& domain = (idom == iRestDom) ? theRestDomElems : theElems[idom];
                   for ( int ivol = 0; ivol < nbvol; ivol++ )
                   {
-                    int smdsId = meshDS->FromVtkToSmds(vtkVolIds[ivol]);
+                    smIdType smdsId = meshDS->FromVtkToSmds(vtkVolIds[ivol]);
                     const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
                     if (domain.count(elem))
                     {
                     const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
                     if (domain.count(elem))
                     {
@@ -11694,7 +11758,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
                       //MESSAGE("  domain " << idom << " volume " << elem->GetID());
                       double values[3] = { 0,0,0 };
                       vtkIdType npts = 0;
                       //MESSAGE("  domain " << idom << " volume " << elem->GetID());
                       double values[3] = { 0,0,0 };
                       vtkIdType npts = 0;
-                      vtkIdType* pts = 0;
+                      vtkIdType const *pts(nullptr);
                       grid->GetCellPoints(vtkVolIds[ivol], npts, pts);
                       for ( vtkIdType i = 0; i < npts; ++i )
                       {
                       grid->GetCellPoints(vtkVolIds[ivol], npts, pts);
                       for ( vtkIdType i = 0; i < npts; ++i )
                       {
@@ -11749,7 +11813,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
 
   //MESSAGE(".. Creation of elements: simple junction");
   std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
 
   //MESSAGE(".. Creation of elements: simple junction");
-  if (createJointElems)
+  if ( createJointElems )
   {
     string joints2DName = "joints2D";
     mapOfJunctionGroups[joints2DName] = this->myMesh->AddGroup(SMDSAbs_Face, joints2DName.c_str());
   {
     string joints2DName = "joints2D";
     mapOfJunctionGroups[joints2DName] = this->myMesh->AddGroup(SMDSAbs_Face, joints2DName.c_str());
@@ -11764,15 +11828,14 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       DownIdType face = itface->first;
       std::set<int> oldNodes;
       std::set<int>::iterator itn;
       DownIdType face = itface->first;
       std::set<int> oldNodes;
       std::set<int>::iterator itn;
-      oldNodes.clear();
       grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
 
       grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
 
-      std::map<int, int> domvol = itface->second;
+      std::map<int, int>          domvol = itface->second;
       std::map<int, int>::iterator itdom = domvol.begin();
       std::map<int, int>::iterator itdom = domvol.begin();
-      int dom1 = itdom->first;
+      int     dom1 = itdom->first;
       int vtkVolId = itdom->second;
       itdom++;
       int vtkVolId = itdom->second;
       itdom++;
-      int dom2 = itdom->first;
+      int           dom2 = itdom->first;
       SMDS_MeshCell *vol = grid->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains,
                                                        nodeQuadDomains);
       stringstream grpname;
       SMDS_MeshCell *vol = grid->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains,
                                                        nodeQuadDomains);
       stringstream grpname;
@@ -11829,7 +11892,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
     std::map<std::vector<int>, std::vector<int> >::iterator ite = edgesMultiDomains.begin();
     for (; ite != edgesMultiDomains.end(); ++ite)
     {
     std::map<std::vector<int>, std::vector<int> >::iterator ite = edgesMultiDomains.begin();
     for (; ite != edgesMultiDomains.end(); ++ite)
     {
-      vector<int> nodes = ite->first;
+      vector<int>    nodes = ite->first;
       vector<int> orderDom = ite->second;
       vector<vtkIdType> orderedNodes;
       if (nodes.size() == 2)
       vector<int> orderDom = ite->second;
       vector<vtkIdType> orderedNodes;
       if (nodes.size() == 2)
@@ -11866,10 +11929,9 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
 
   std::map<DownIdType, std::map<int,int>, DownIdCompare> faceOrEdgeDom; // cellToModify --> (id domain --> id cell)
   std::map<int,int> feDom; // vtk id of cell to modify --> id domain
 
   std::map<DownIdType, std::map<int,int>, DownIdCompare> faceOrEdgeDom; // cellToModify --> (id domain --> id cell)
   std::map<int,int> feDom; // vtk id of cell to modify --> id domain
-  faceOrEdgeDom.clear();
-  feDom.clear();
 
   //MESSAGE(".. Modification of elements");
 
   //MESSAGE(".. Modification of elements");
+  SMDSAbs_ElementType domainType = (*theElems[0].begin())->GetType();
   for (int idomain = idom0; idomain < nbDomains; idomain++)
   {
     std::map<int, std::map<int, int> >::const_iterator itnod = nodeDomains.begin();
   for (int idomain = idom0; idomain < nbDomains; idomain++)
   {
     std::map<int, std::map<int, int> >::const_iterator itnod = nodeDomains.begin();
@@ -11877,7 +11939,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
     {
       int oldId = itnod->first;
       //MESSAGE("     node " << oldId);
     {
       int oldId = itnod->first;
       //MESSAGE("     node " << oldId);
-      vtkCellLinks::Link l = grid->GetCellLinks()->GetLink(oldId);
+      vtkCellLinks::Link l = (static_cast< vtkCellLinks *>(grid->GetCellLinks()))->GetLink(oldId);
       for (int i = 0; i < l.ncells; i++)
       {
         int vtkId = l.cells[i];
       for (int i = 0; i < l.ncells; i++)
       {
         int vtkId = l.cells[i];
@@ -11887,13 +11949,29 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
           continue; // new cells: not to be modified
         DownIdType aCell(downId, vtkType);
         int volParents[1000];
           continue; // new cells: not to be modified
         DownIdType aCell(downId, vtkType);
         int volParents[1000];
-        int nbvol = grid->GetParentVolumes(volParents, vtkId);
+        int nbvol = 0;
+        nbvol = grid->GetParentVolumes(volParents, vtkId);
+        if ( domainType == SMDSAbs_Volume )
+        {
+          nbvol = grid->GetParentVolumes(volParents, vtkId);
+        }
+        else // domainType == SMDSAbs_Face
+        {
+          const int            nbFaces = grid->getDownArray(vtkType)->getNumberOfUpCells(downId);
+          const int           *upCells = grid->getDownArray(vtkType)->getUpCells(downId);
+          const unsigned char* upTypes = grid->getDownArray(vtkType)->getUpTypes(downId);
+          for (int i=0; i< nbFaces; i++)
+          {
+            int vtkFaceId = grid->getDownArray( upTypes[i] )->getVtkCellId(upCells[i]);
+            if (vtkFaceId >= 0)
+              volParents[nbvol++] = vtkFaceId;
+          }
+        }
         for (int j = 0; j < nbvol; j++)
           if (celldom.count(volParents[j]) && (celldom[volParents[j]] == idomain))
             if (!feDom.count(vtkId))
             {
               feDom[vtkId] = idomain;
         for (int j = 0; j < nbvol; j++)
           if (celldom.count(volParents[j]) && (celldom[volParents[j]] == idomain))
             if (!feDom.count(vtkId))
             {
               feDom[vtkId] = idomain;
-              faceOrEdgeDom[aCell] = emptyMap;
               faceOrEdgeDom[aCell][idomain] = vtkId; // affect face or edge to the first domain only
               //MESSAGE("affect cell " << this->GetMeshDS()->FromVtkToSmds(vtkId) << " domain " << idomain
               //        << " type " << vtkType << " downId " << downId);
               faceOrEdgeDom[aCell][idomain] = vtkId; // affect face or edge to the first domain only
               //MESSAGE("affect cell " << this->GetMeshDS()->FromVtkToSmds(vtkId) << " domain " << idomain
               //        << " type " << vtkType << " downId " << downId);
@@ -11916,7 +11994,6 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       DownIdType face = itface->first;
       std::set<int> oldNodes;
       std::set<int>::iterator itn;
       DownIdType face = itface->first;
       std::set<int> oldNodes;
       std::set<int>::iterator itn;
-      oldNodes.clear();
       grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
       //MESSAGE("examine cell, downId " << face.cellId << " type " << int(face.cellType));
       std::map<int, int> localClonedNodeIds;
       grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
       //MESSAGE("examine cell, downId " << face.cellId << " type " << int(face.cellType));
       std::map<int, int> localClonedNodeIds;
@@ -11983,10 +12060,7 @@ bool SMESH_MeshEditor::CreateFlatElementsOnFacesGroups(const std::vector<TIDSort
 
   std::map<const SMDS_MeshNode*, const SMDS_MeshNode*> clonedNodes;
   std::map<const SMDS_MeshNode*, const SMDS_MeshNode*> intermediateNodes;
 
   std::map<const SMDS_MeshNode*, const SMDS_MeshNode*> clonedNodes;
   std::map<const SMDS_MeshNode*, const SMDS_MeshNode*> intermediateNodes;
-  clonedNodes.clear();
-  intermediateNodes.clear();
   std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
   std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
-  mapOfJunctionGroups.clear();
 
   for ( size_t idom = 0; idom < theElems.size(); idom++ )
   {
 
   for ( size_t idom = 0; idom < theElems.size(); idom++ )
   {
@@ -12233,7 +12307,6 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
   std::set<int> setOfVolToCheck;
 
   std::vector<gp_Pnt> gpnts;
   std::set<int> setOfVolToCheck;
 
   std::vector<gp_Pnt> gpnts;
-  gpnts.clear();
 
   if (isNodeGroup) // --- a group of nodes is provided : find all the volumes using one or more of this nodes
   {
 
   if (isNodeGroup) // --- a group of nodes is provided : find all the volumes using one or more of this nodes
   {
@@ -12320,7 +12393,6 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
     //     Fill the group of inside volumes
 
     std::map<int, double> mapOfNodeDistance2;
     //     Fill the group of inside volumes
 
     std::map<int, double> mapOfNodeDistance2;
-    mapOfNodeDistance2.clear();
     std::set<int> setOfOutsideVol;
     while (!setOfVolToCheck.empty())
     {
     std::set<int> setOfOutsideVol;
     while (!setOfVolToCheck.empty())
     {
@@ -12329,7 +12401,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
       //MESSAGE("volume to check,  vtkId " << vtkId << " smdsId " << meshDS->FromVtkToSmds(vtkId));
       bool volInside = false;
       vtkIdType npts = 0;
       //MESSAGE("volume to check,  vtkId " << vtkId << " smdsId " << meshDS->FromVtkToSmds(vtkId));
       bool volInside = false;
       vtkIdType npts = 0;
-      vtkIdType* pts = 0;
+      vtkIdType const *pts(nullptr);
       grid->GetCellPoints(vtkId, npts, pts);
       for (int i=0; i<npts; i++)
       {
       grid->GetCellPoints(vtkId, npts, pts);
       for (int i=0; i<npts; i++)
       {
@@ -12497,9 +12569,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
   //     project polylines on subshapes, and partition, to get geom faces
 
   std::map<int, std::set<int> > shapeIdToVtkIdSet; // shapeId --> set of vtkId on skin
   //     project polylines on subshapes, and partition, to get geom faces
 
   std::map<int, std::set<int> > shapeIdToVtkIdSet; // shapeId --> set of vtkId on skin
-  std::set<int> emptySet;
-  emptySet.clear();
-  std::set<int> shapeIds;
+  std::set<int>                 shapeIds;
 
   SMDS_ElemIteratorPtr itelem = sgrps->GetElements();
   while (itelem->more())
 
   SMDS_ElemIteratorPtr itelem = sgrps->GetElements();
   while (itelem->more())
@@ -12509,7 +12579,6 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
     int   vtkId = elem->GetVtkID();
     if (!shapeIdToVtkIdSet.count(shapeId))
     {
     int   vtkId = elem->GetVtkID();
     if (!shapeIdToVtkIdSet.count(shapeId))
     {
-      shapeIdToVtkIdSet[shapeId] = emptySet;
       shapeIds.insert(shapeId);
     }
     shapeIdToVtkIdSet[shapeId].insert(vtkId);
       shapeIds.insert(shapeId);
     }
     shapeIdToVtkIdSet[shapeId].insert(vtkId);
@@ -12517,7 +12586,6 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
 
   std::map<int, std::set<DownIdType, DownIdCompare> > shapeIdToEdges; // shapeId --> set of downward edges
   std::set<DownIdType, DownIdCompare> emptyEdges;
 
   std::map<int, std::set<DownIdType, DownIdCompare> > shapeIdToEdges; // shapeId --> set of downward edges
   std::set<DownIdType, DownIdCompare> emptyEdges;
-  emptyEdges.clear();
 
   std::map<int, std::set<int> >::iterator itShape =  shapeIdToVtkIdSet.begin();
   for (; itShape != shapeIdToVtkIdSet.end(); ++itShape)
 
   std::map<int, std::set<int> >::iterator itShape =  shapeIdToVtkIdSet.begin();
   for (; itShape != shapeIdToVtkIdSet.end(); ++itShape)
@@ -12541,7 +12609,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
       {
         if (neighborsVtkIds[n]<0) // only smds faces are considered as neighbors here
           continue;
       {
         if (neighborsVtkIds[n]<0) // only smds faces are considered as neighbors here
           continue;
-        int smdsId = meshDS->FromVtkToSmds(neighborsVtkIds[n]);
+        smIdType smdsId = meshDS->FromVtkToSmds(neighborsVtkIds[n]);
         const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
         if ( shapeIds.count(elem->getshapeId()) && !sgrps->Contains(elem)) // edge : neighbor in the set of shape, not in the group
         {
         const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
         if ( shapeIds.count(elem->getshapeId()) && !sgrps->Contains(elem)) // edge : neighbor in the set of shape, not in the group
         {
@@ -12560,7 +12628,6 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
     }
 
     std::list<int> order;
     }
 
     std::list<int> order;
-    order.clear();
     if (nodesEdges.size() > 0)
     {
       order.push_back(nodesEdges[0]); //MESSAGE("       --- back " << order.back()+1); // SMDS id = VTK id + 1;
     if (nodesEdges.size() > 0)
     {
       order.push_back(nodesEdges[0]); //MESSAGE("       --- back " << order.back()+1); // SMDS id = VTK id + 1;
@@ -12896,8 +12963,8 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
         for ( inode = 0; inode < srcNodes.size(); ++inode )
           tgtNodes[inode] = getNodeWithSameID( tgtMeshDS, srcNodes[inode] );
         if ( /*aroundElements && */tgtEditor.GetMeshDS()->FindElement( tgtNodes,
         for ( inode = 0; inode < srcNodes.size(); ++inode )
           tgtNodes[inode] = getNodeWithSameID( tgtMeshDS, srcNodes[inode] );
         if ( /*aroundElements && */tgtEditor.GetMeshDS()->FindElement( tgtNodes,
-                                                                   missType,
-                                                                   /*noMedium=*/false))
+                                                                       missType,
+                                                                       /*noMedium=*/false))
           continue;
         tgtEditor.AddElement( tgtNodes, elemKind.SetPoly( tgtNodes.size()/(iQuad+1) > 4 ));
         ++nbAddedBnd;
           continue;
         tgtEditor.AddElement( tgtNodes, elemKind.SetPoly( tgtNodes.size()/(iQuad+1) > 4 ));
         ++nbAddedBnd;
@@ -12907,8 +12974,8 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
       {
         TConnectivity& nodes = missingBndElems[ i ];
         if ( /*aroundElements && */tgtEditor.GetMeshDS()->FindElement( nodes,
       {
         TConnectivity& nodes = missingBndElems[ i ];
         if ( /*aroundElements && */tgtEditor.GetMeshDS()->FindElement( nodes,
-                                                                   missType,
-                                                                   /*noMedium=*/false))
+                                                                       missType,
+                                                                       /*noMedium=*/false))
           continue;
         SMDS_MeshElement* newElem =
           tgtEditor.AddElement( nodes, elemKind.SetPoly( nodes.size()/(iQuad+1) > 4 ));
           continue;
         SMDS_MeshElement* newElem =
           tgtEditor.AddElement( nodes, elemKind.SetPoly( nodes.size()/(iQuad+1) > 4 ));