Salome HOME
23050: EDF 10631 SMESH: Nodes outside the plane after ConvertToQuadratic
[modules/smesh.git] / src / SMESH / SMESH_MeshEditor.cxx
index 6bdd80cac2d2dcd5c032006599517cfb884ebeba..436a0ad0339816e6c832f0ef080d07123783653c 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2014  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2015  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
@@ -50,6 +50,7 @@
 #include <Basics_OCCTVersion.hxx>
 
 #include "utilities.h"
+#include "chrono.hxx"
 
 #include <BRepAdaptor_Surface.hxx>
 #include <BRepBuilderAPI_MakeEdge.hxx>
@@ -130,7 +131,7 @@ SMESH_MeshEditor::SMESH_MeshEditor( SMESH_Mesh* theMesh )
  */
 //================================================================================
 
-void SMESH_MeshEditor::CrearLastCreated()
+void SMESH_MeshEditor::ClearLastCreated()
 {
   myLastCreatedNodes.Clear();
   myLastCreatedElems.Clear();
@@ -1243,6 +1244,92 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
   return nbReori;
 }
 
+//================================================================================
+/*!
+ * \brief Reorient faces basing on orientation of adjacent volumes.
+ * \param theFaces - faces to reorient. If empty, all mesh faces are treated.
+ * \param theVolumes - reference volumes.
+ * \param theOutsideNormal - to orient faces to have their normal
+ *        pointing either \a outside or \a inside the adjacent volumes.
+ * \return number of reoriented faces.
+ */
+//================================================================================
+
+int SMESH_MeshEditor::Reorient2DBy3D (TIDSortedElemSet & theFaces,
+                                      TIDSortedElemSet & theVolumes,
+                                      const bool         theOutsideNormal)
+{
+  int nbReori = 0;
+
+  SMDS_ElemIteratorPtr faceIt;
+  if ( theFaces.empty() )
+    faceIt = GetMeshDS()->elementsIterator( SMDSAbs_Face );
+  else
+    faceIt = elemSetIterator( theFaces );
+
+  vector< const SMDS_MeshNode* > faceNodes;
+  TIDSortedElemSet checkedVolumes;
+  set< const SMDS_MeshNode* > faceNodesSet;
+  SMDS_VolumeTool volumeTool;
+
+  while ( faceIt->more() ) // loop on given faces
+  {
+    const SMDS_MeshElement* face = faceIt->next();
+    if ( face->GetType() != SMDSAbs_Face )
+      continue;
+
+    const int nbCornersNodes = face->NbCornerNodes();
+    faceNodes.assign( face->begin_nodes(), face->end_nodes() );
+
+    checkedVolumes.clear();
+    SMDS_ElemIteratorPtr vIt = faceNodes[ 0 ]->GetInverseElementIterator( SMDSAbs_Volume );
+    while ( vIt->more() )
+    {
+      const SMDS_MeshElement* volume = vIt->next();
+
+      if ( !checkedVolumes.insert( volume ).second )
+        continue;
+      if ( !theVolumes.empty() && !theVolumes.count( volume ))
+        continue;
+
+      // is volume adjacent?
+      bool allNodesCommon = true;
+      for ( int iN = 1; iN < nbCornersNodes && allNodesCommon; ++iN )
+        allNodesCommon = ( volume->GetNodeIndex( faceNodes[ iN ]) > -1 );
+      if ( !allNodesCommon )
+        continue;
+
+      // get nodes of a corresponding volume facet
+      faceNodesSet.clear();
+      faceNodesSet.insert( faceNodes.begin(), faceNodes.end() );
+      volumeTool.Set( volume );
+      int facetID = volumeTool.GetFaceIndex( faceNodesSet );
+      if ( facetID < 0 ) continue;
+      volumeTool.SetExternalNormal();
+      const SMDS_MeshNode** facetNodes = volumeTool.GetFaceNodes( facetID );
+
+      // compare order of faceNodes and facetNodes
+      const int iQ = 1 + ( nbCornersNodes < faceNodes.size() );
+      int iNN[2];
+      for ( int i = 0; i < 2; ++i )
+      {
+        const SMDS_MeshNode* n = facetNodes[ i*iQ ];
+        for ( int iN = 0; iN < nbCornersNodes; ++iN )
+          if ( faceNodes[ iN ] == n )
+          {
+            iNN[ i ] = iN;
+            break;
+          }
+      }
+      bool isOutside = Abs( iNN[0]-iNN[1] ) == 1 ? iNN[0] < iNN[1] : iNN[0] > iNN[1];
+      if ( isOutside != theOutsideNormal )
+        nbReori += Reorient( face );
+    }
+  }  // loop on given faces
+
+  return nbReori;
+}
+
 //=======================================================================
 //function : getBadRate
 //purpose  :
@@ -2093,10 +2180,6 @@ namespace
 void SMESH_MeshEditor::SplitVolumes (const TFacetOfElem & theElems,
                                      const int            theMethodFlags)
 {
-  // std-like iterator on coordinates of nodes of mesh element
-  typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > NXyzIterator;
-  NXyzIterator xyzEnd;
-
   SMDS_VolumeTool    volTool;
   SMESH_MesherHelper helper( *GetMesh()), fHelper(*GetMesh());
   fHelper.ToFixNodeParameters( true );
@@ -3615,13 +3698,8 @@ static bool getClosestUV (Extrema_GenExtPS& projector,
   if ( projector.IsDone() ) {
     double u, v, minVal = DBL_MAX;
     for ( int i = projector.NbExt(); i > 0; i-- )
-#if OCC_VERSION_LARGE > 0x06040000 // Porting to OCCT6.5.1
       if ( projector.SquareDistance( i ) < minVal ) {
         minVal = projector.SquareDistance( i );
-#else
-      if ( projector.Value( i ) < minVal ) {
-        minVal = projector.Value( i );
-#endif
         projector.Point( i ).Parameter( u, v );
       }
     result.SetCoord( u, v );
@@ -4138,26 +4216,49 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
 
 }
 
-//=======================================================================
-//function : isReverse
-//purpose  : Return true if normal of prevNodes is not co-directied with
-//           gp_Vec(prevNodes[iNotSame],nextNodes[iNotSame]).
-//           iNotSame is where prevNodes and nextNodes are different.
-//           If result is true then future volume orientation is OK
-//=======================================================================
-
-static bool isReverse(const SMDS_MeshElement*             face,
-                      const vector<const SMDS_MeshNode*>& prevNodes,
-                      const vector<const SMDS_MeshNode*>& nextNodes,
-                      const int                           iNotSame)
+namespace
 {
+  //=======================================================================
+  //function : isReverse
+  //purpose  : Return true if normal of prevNodes is not co-directied with
+  //           gp_Vec(prevNodes[iNotSame],nextNodes[iNotSame]).
+  //           iNotSame is where prevNodes and nextNodes are different.
+  //           If result is true then future volume orientation is OK
+  //=======================================================================
 
-  SMESH_TNodeXYZ pP = prevNodes[ iNotSame ];
-  SMESH_TNodeXYZ pN = nextNodes[ iNotSame ];
-  gp_XYZ extrDir( pN - pP ), faceNorm;
-  SMESH_MeshAlgos::FaceNormal( face, faceNorm, /*normalized=*/false );
+  bool isReverse(const SMDS_MeshElement*             face,
+                 const vector<const SMDS_MeshNode*>& prevNodes,
+                 const vector<const SMDS_MeshNode*>& nextNodes,
+                 const int                           iNotSame)
+  {
+
+    SMESH_TNodeXYZ pP = prevNodes[ iNotSame ];
+    SMESH_TNodeXYZ pN = nextNodes[ iNotSame ];
+    gp_XYZ extrDir( pN - pP ), faceNorm;
+    SMESH_MeshAlgos::FaceNormal( face, faceNorm, /*normalized=*/false );
+
+    return faceNorm * extrDir < 0.0;
+  }
 
-  return faceNorm * extrDir < 0.0;
+  //================================================================================
+  /*!
+   * \brief Assure that theElemSets[0] holds elements, not nodes
+   */
+  //================================================================================
+
+  void setElemsFirst( TIDSortedElemSet theElemSets[2] )
+  {
+    if ( !theElemSets[0].empty() &&
+         (*theElemSets[0].begin())->GetType() == SMDSAbs_Node )
+    {
+      std::swap( theElemSets[0], theElemSets[1] );
+    }
+    else if ( !theElemSets[1].empty() &&
+              (*theElemSets[1].begin())->GetType() != SMDSAbs_Node )
+    {
+      std::swap( theElemSets[0], theElemSets[1] );
+    }
+  }
 }
 
 //=======================================================================
@@ -4257,6 +4358,30 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
       }
     }
   }
+  else if ( elem->GetType() == SMDSAbs_Edge )
+  {
+    // orient a new face same as adjacent one
+    int i1, i2;
+    const SMDS_MeshElement* e;
+    TIDSortedElemSet dummy;
+    if (( e = SMESH_MeshAlgos::FindFaceInSet( nextNod[0], prevNod[0], dummy,dummy, &i1, &i2 )) ||
+        ( e = SMESH_MeshAlgos::FindFaceInSet( prevNod[1], nextNod[1], dummy,dummy, &i1, &i2 )) ||
+        ( e = SMESH_MeshAlgos::FindFaceInSet( prevNod[0], prevNod[1], dummy,dummy, &i1, &i2 )))
+    {
+      // there is an adjacent face, check order of nodes in it
+      bool sameOrder = ( Abs( i2 - i1 ) == 1 ) ? ( i2 > i1 ) : ( i2 < i1 );
+      if ( sameOrder )
+      {
+        std::swap( itNN[0],    itNN[1] );
+        std::swap( prevNod[0], prevNod[1] );
+        std::swap( nextNod[0], nextNod[1] );
+        isSingleNode.swap( isSingleNode[0], isSingleNode[1] );
+        if ( nbSame > 0 )
+          sames[0] = 1 - sames[0];
+        iNotSameNode = 1 - iNotSameNode;
+      }
+    }
+  }
 
   int iSameNode = 0, iBeforeSame = 0, iAfterSame = 0, iOpposSame = 0;
   if ( nbSame > 0 ) {
@@ -4362,11 +4487,11 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
                 return; // medium node on axis
               }
               else if(sames[0]==0)
-                aNewElem = aMesh->AddFace(prevNod[0], nextNod[1], prevNod[1],
-                                          nextNod[2], midlNod[1], prevNod[2]);
+                aNewElem = aMesh->AddFace(prevNod[0], prevNod[1], nextNod[1],
+                                          prevNod[2], midlNod[1], nextNod[2] );
               else // sames[0]==1
-                aNewElem = aMesh->AddFace(prevNod[0], nextNod[0], prevNod[1],
-                                          midlNod[0], nextNod[2], prevNod[2]);
+                aNewElem = aMesh->AddFace(prevNod[0], prevNod[1], nextNod[0],
+                                          prevNod[2], nextNod[2], midlNod[0]);
             }
           }
           else if ( nbDouble == 3 )
@@ -4532,8 +4657,8 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
 
       default:
         break;
-      }
-    }
+      } // switch ( baseType )
+    } // scope
 
     if ( !aNewElem && elem->GetType() == SMDSAbs_Face ) // try to create a polyherdal prism
     {
@@ -4587,7 +4712,8 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
           polyedre_nodes.resize( prevNbNodes );
       }
       aNewElem = aMesh->AddPolyhedralVolume (polyedre_nodes, quantities);
-    }
+
+    } //  // try to create a polyherdal prism
 
     if ( aNewElem ) {
       newElems.push_back( aNewElem );
@@ -4599,7 +4725,7 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
     for ( iNode = 0; iNode < nbNodes; iNode++ )
       prevNod[ iNode ] = nextNod[ iNode ];
 
-  } // for steps
+  } // loop on steps
 }
 
 //=======================================================================
@@ -5018,7 +5144,7 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
 //=======================================================================
 
 SMESH_MeshEditor::PGroupIDs
-SMESH_MeshEditor::RotationSweep(TIDSortedElemSet & theElems,
+SMESH_MeshEditor::RotationSweep(TIDSortedElemSet   theElemSets[2],
                                 const gp_Ax1&      theAxis,
                                 const double       theAngle,
                                 const int          theNbSteps,
@@ -5050,84 +5176,89 @@ SMESH_MeshEditor::RotationSweep(TIDSortedElemSet & theElems,
   const bool isQuadraticMesh = bool( myMesh->NbEdges(ORDER_QUADRATIC) +
                                      myMesh->NbFaces(ORDER_QUADRATIC) +
                                      myMesh->NbVolumes(ORDER_QUADRATIC) );
-  // loop on theElems
+  // loop on theElemSets
+  setElemsFirst( theElemSets );
   TIDSortedElemSet::iterator itElem;
-  for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
-    const SMDS_MeshElement* elem = *itElem;
-    if ( !elem || elem->GetType() == SMDSAbs_Volume )
-      continue;
-    vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
-    newNodesItVec.reserve( elem->NbNodes() );
+  for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
+  {
+    TIDSortedElemSet& theElems = theElemSets[ is2ndSet ];
+    for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
+      const SMDS_MeshElement* elem = *itElem;
+      if ( !elem || elem->GetType() == SMDSAbs_Volume )
+        continue;
+      vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
+      newNodesItVec.reserve( elem->NbNodes() );
 
-    // loop on elem nodes
-    SMDS_ElemIteratorPtr itN = elem->nodesIterator();
-    while ( itN->more() )
-    {
-      // check if a node has been already sweeped
-      const SMDS_MeshNode* node = cast2Node( itN->next() );
+      // loop on elem nodes
+      SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+      while ( itN->more() )
+      {
+        // check if a node has been already sweeped
+        const SMDS_MeshNode* node = cast2Node( itN->next() );
 
-      gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
-      double coord[3];
-      aXYZ.Coord( coord[0], coord[1], coord[2] );
-      bool isOnAxis = ( aLine.SquareDistance( aXYZ ) <= aSqTol );
+        gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
+        double coord[3];
+        aXYZ.Coord( coord[0], coord[1], coord[2] );
+        bool isOnAxis = ( aLine.SquareDistance( aXYZ ) <= aSqTol );
 
-      TNodeOfNodeListMapItr nIt =
-        mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
-      list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
-      if ( listNewNodes.empty() )
-      {
-        // check if we are to create medium nodes between corner ones
-        bool needMediumNodes = false;
-        if ( isQuadraticMesh )
+        TNodeOfNodeListMapItr nIt =
+          mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
+        list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
+        if ( listNewNodes.empty() )
         {
-          SMDS_ElemIteratorPtr it = node->GetInverseElementIterator();
-          while (it->more() && !needMediumNodes )
+          // check if we are to create medium nodes between corner ones
+          bool needMediumNodes = false;
+          if ( isQuadraticMesh )
           {
-            const SMDS_MeshElement* invElem = it->next();
-            if ( invElem != elem && !theElems.count( invElem )) continue;
-            needMediumNodes = ( invElem->IsQuadratic() && !invElem->IsMediumNode(node) );
-            if ( !needMediumNodes && invElem->GetEntityType() == SMDSEntity_BiQuad_Quadrangle )
-              needMediumNodes = true;
+            SMDS_ElemIteratorPtr it = node->GetInverseElementIterator();
+            while (it->more() && !needMediumNodes )
+            {
+              const SMDS_MeshElement* invElem = it->next();
+              if ( invElem != elem && !theElems.count( invElem )) continue;
+              needMediumNodes = ( invElem->IsQuadratic() && !invElem->IsMediumNode(node) );
+              if ( !needMediumNodes && invElem->GetEntityType() == SMDSEntity_BiQuad_Quadrangle )
+                needMediumNodes = true;
+            }
           }
-        }
 
-        // make new nodes
-        const SMDS_MeshNode * newNode = node;
-        for ( int i = 0; i < theNbSteps; i++ ) {
-          if ( !isOnAxis ) {
-            if ( needMediumNodes )  // create a medium node
-            {
-              aTrsf2.Transforms( coord[0], coord[1], coord[2] );
+          // make new nodes
+          const SMDS_MeshNode * newNode = node;
+          for ( int i = 0; i < theNbSteps; i++ ) {
+            if ( !isOnAxis ) {
+              if ( needMediumNodes )  // create a medium node
+              {
+                aTrsf2.Transforms( coord[0], coord[1], coord[2] );
+                newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
+                myLastCreatedNodes.Append(newNode);
+                srcNodes.Append( node );
+                listNewNodes.push_back( newNode );
+                aTrsf2.Transforms( coord[0], coord[1], coord[2] );
+              }
+              else {
+                aTrsf.Transforms( coord[0], coord[1], coord[2] );
+              }
+              // create a corner node
               newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
               myLastCreatedNodes.Append(newNode);
               srcNodes.Append( node );
               listNewNodes.push_back( newNode );
-              aTrsf2.Transforms( coord[0], coord[1], coord[2] );
             }
             else {
-              aTrsf.Transforms( coord[0], coord[1], coord[2] );
+              listNewNodes.push_back( newNode );
+              // if ( needMediumNodes )
+              //   listNewNodes.push_back( newNode );
             }
-            // create a corner node
-            newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
-            myLastCreatedNodes.Append(newNode);
-            srcNodes.Append( node );
-            listNewNodes.push_back( newNode );
-          }
-          else {
-            listNewNodes.push_back( newNode );
-            // if ( needMediumNodes )
-            //   listNewNodes.push_back( newNode );
           }
         }
+        newNodesItVec.push_back( nIt );
       }
-      newNodesItVec.push_back( nIt );
+      // make new elements
+      sweepElement( elem, newNodesItVec, newElemsMap[elem], theNbSteps, srcElems );
     }
-    // make new elements
-    sweepElement( elem, newNodesItVec, newElemsMap[elem], theNbSteps, srcElems );
   }
 
   if ( theMakeWalls )
-    makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElems, theNbSteps, srcElems );
+    makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElemSets[0], theNbSteps, srcElems );
 
   PGroupIDs newGroupIDs;
   if ( theMakeGroups )
@@ -5136,50 +5267,328 @@ SMESH_MeshEditor::RotationSweep(TIDSortedElemSet & theElems,
   return newGroupIDs;
 }
 
+//=======================================================================
+//function : ExtrusParam
+//purpose  : standard construction
+//=======================================================================
+
+SMESH_MeshEditor::ExtrusParam::ExtrusParam( const gp_Vec&  theStep,
+                                            const int      theNbSteps,
+                                            const int      theFlags,
+                                            const double   theTolerance):
+  myDir( theStep ),
+  myFlags( theFlags ),
+  myTolerance( theTolerance ),
+  myElemsToUse( NULL )
+{
+  mySteps = new TColStd_HSequenceOfReal;
+  const double stepSize = theStep.Magnitude();
+  for (int i=1; i<=theNbSteps; i++ )
+    mySteps->Append( stepSize );
+
+  if (( theFlags & EXTRUSION_FLAG_SEW ) &&
+      ( theTolerance > 0 ))
+  {
+    myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByDirAndSew;
+  }
+  else
+  {
+    myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByDir;
+  }
+}
 
 //=======================================================================
-//function : CreateNode
-//purpose  :
+//function : ExtrusParam
+//purpose  : steps are given explicitly
 //=======================================================================
-const SMDS_MeshNode* SMESH_MeshEditor::CreateNode(const double x,
-                                                  const double y,
-                                                  const double z,
-                                                  const double tolnode,
-                                                  SMESH_SequenceOfNode& aNodes)
+
+SMESH_MeshEditor::ExtrusParam::ExtrusParam( const gp_Dir&                   theDir,
+                                            Handle(TColStd_HSequenceOfReal) theSteps,
+                                            const int                       theFlags,
+                                            const double                    theTolerance):
+  myDir( theDir ),
+  mySteps( theSteps ),
+  myFlags( theFlags ),
+  myTolerance( theTolerance ),
+  myElemsToUse( NULL )
 {
-  // myLastCreatedElems.Clear();
-  // myLastCreatedNodes.Clear();
+  if (( theFlags & EXTRUSION_FLAG_SEW ) &&
+      ( theTolerance > 0 ))
+  {
+    myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByDirAndSew;
+  }
+  else
+  {
+    myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByDir;
+  }
+}
 
-  gp_Pnt P1(x,y,z);
-  SMESHDS_Mesh * aMesh = myMesh->GetMeshDS();
+//=======================================================================
+//function : ExtrusParam
+//purpose  : for extrusion by normal
+//=======================================================================
+
+SMESH_MeshEditor::ExtrusParam::ExtrusParam( const double theStepSize,
+                                            const int    theNbSteps,
+                                            const int    theFlags,
+                                            const int    theDim ):
+  myDir( 1,0,0 ),
+  mySteps( new TColStd_HSequenceOfReal ),
+  myFlags( theFlags ),
+  myTolerance( 0 ),
+  myElemsToUse( NULL )
+{
+  for (int i = 0; i < theNbSteps; i++ )
+    mySteps->Append( theStepSize );
+
+  if ( theDim == 1 )
+  {
+    myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByNormal1D;
+  }
+  else
+  {
+    myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByNormal2D;
+  }
+}
 
-  // try to search in sequence of existing nodes
-  // if aNodes.Length()>0 we 'nave to use given sequence
-  // else - use all nodes of mesh
-  if(aNodes.Length()>0) {
-    int i;
-    for(i=1; i<=aNodes.Length(); i++) {
-      gp_Pnt P2(aNodes.Value(i)->X(),aNodes.Value(i)->Y(),aNodes.Value(i)->Z());
-      if(P1.Distance(P2)<tolnode)
-        return aNodes.Value(i);
+//=======================================================================
+//function : ExtrusParam::SetElementsToUse
+//purpose  : stores elements to use for extrusion by normal, depending on
+//           state of EXTRUSION_FLAG_USE_INPUT_ELEMS_ONLY flag
+//=======================================================================
+
+void SMESH_MeshEditor::ExtrusParam::SetElementsToUse( const TIDSortedElemSet& elems )
+{
+  myElemsToUse = ToUseInpElemsOnly() ? & elems : 0;
+}
+
+//=======================================================================
+//function : ExtrusParam::beginStepIter
+//purpose  : prepare iteration on steps
+//=======================================================================
+
+void SMESH_MeshEditor::ExtrusParam::beginStepIter( bool withMediumNodes )
+{
+  myWithMediumNodes = withMediumNodes;
+  myNextStep = 1;
+  myCurSteps.clear();
+}
+//=======================================================================
+//function : ExtrusParam::moreSteps
+//purpose  : are there more steps?
+//=======================================================================
+
+bool SMESH_MeshEditor::ExtrusParam::moreSteps()
+{
+  return myNextStep <= mySteps->Length() || !myCurSteps.empty();
+}
+//=======================================================================
+//function : ExtrusParam::nextStep
+//purpose  : returns the next step
+//=======================================================================
+
+double SMESH_MeshEditor::ExtrusParam::nextStep()
+{
+  double res = 0;
+  if ( !myCurSteps.empty() )
+  {
+    res = myCurSteps.back();
+    myCurSteps.pop_back();
+  }
+  else if ( myNextStep <= mySteps->Length() )
+  {
+    myCurSteps.push_back( mySteps->Value( myNextStep ));
+    ++myNextStep;
+    if ( myWithMediumNodes )
+    {
+      myCurSteps.back() /= 2.;
+      myCurSteps.push_back( myCurSteps.back() );
     }
+    res = nextStep();
   }
-  else {
-    SMDS_NodeIteratorPtr itn = aMesh->nodesIterator();
-    while(itn->more()) {
-      const SMDS_MeshNode* aN = static_cast<const SMDS_MeshNode*> (itn->next());
-      gp_Pnt P2(aN->X(),aN->Y(),aN->Z());
-      if(P1.Distance(P2)<tolnode)
-        return aN;
+  return res;
+}
+
+//=======================================================================
+//function : ExtrusParam::makeNodesByDir
+//purpose  : create nodes for standard extrusion
+//=======================================================================
+
+int SMESH_MeshEditor::ExtrusParam::
+makeNodesByDir( SMESHDS_Mesh*                     mesh,
+                const SMDS_MeshNode*              srcNode,
+                std::list<const SMDS_MeshNode*> & newNodes,
+                const bool                        makeMediumNodes)
+{
+  gp_XYZ p = SMESH_TNodeXYZ( srcNode );
+
+  int nbNodes = 0;
+  for ( beginStepIter( makeMediumNodes ); moreSteps(); ++nbNodes ) // loop on steps
+  {
+    p += myDir.XYZ() * nextStep();
+    const SMDS_MeshNode * newNode = mesh->AddNode( p.X(), p.Y(), p.Z() );
+    newNodes.push_back( newNode );
+  }
+  return nbNodes;
+}
+
+//=======================================================================
+//function : ExtrusParam::makeNodesByDirAndSew
+//purpose  : create nodes for standard extrusion with sewing
+//=======================================================================
+
+int SMESH_MeshEditor::ExtrusParam::
+makeNodesByDirAndSew( SMESHDS_Mesh*                     mesh,
+                      const SMDS_MeshNode*              srcNode,
+                      std::list<const SMDS_MeshNode*> & newNodes,
+                      const bool                        makeMediumNodes)
+{
+  gp_XYZ P1 = SMESH_TNodeXYZ( srcNode );
+
+  int nbNodes = 0;
+  for ( beginStepIter( makeMediumNodes ); moreSteps(); ++nbNodes ) // loop on steps
+  {
+    P1 += myDir.XYZ() * nextStep();
+
+    // try to search in sequence of existing nodes
+    // if myNodes.Length()>0 we 'nave to use given sequence
+    // else - use all nodes of mesh
+    const SMDS_MeshNode * node = 0;
+    if ( myNodes.Length() > 0 ) {
+      int i;
+      for(i=1; i<=myNodes.Length(); i++) {
+        gp_XYZ P2 = SMESH_TNodeXYZ( myNodes.Value(i) );
+        if (( P1 - P2 ).SquareModulus() < myTolerance * myTolerance )
+        {
+          node = myNodes.Value(i);
+          break;
+        }
+      }
+    }
+    else {
+      SMDS_NodeIteratorPtr itn = mesh->nodesIterator();
+      while(itn->more()) {
+        SMESH_TNodeXYZ P2( itn->next() );
+        if (( P1 - P2 ).SquareModulus() < myTolerance * myTolerance )
+        {
+          node = P2._node;
+          break;
+        }
+      }
+    }
+
+    if ( !node )
+      node = mesh->AddNode( P1.X(), P1.Y(), P1.Z() );
+
+    newNodes.push_back( node );
+
+  } // loop on steps
+
+  return nbNodes;
+}
+
+//=======================================================================
+//function : ExtrusParam::makeNodesByNormal2D
+//purpose  : create nodes for extrusion using normals of faces
+//=======================================================================
+
+int SMESH_MeshEditor::ExtrusParam::
+makeNodesByNormal2D( SMESHDS_Mesh*                     mesh,
+                     const SMDS_MeshNode*              srcNode,
+                     std::list<const SMDS_MeshNode*> & newNodes,
+                     const bool                        makeMediumNodes)
+{
+  const bool alongAvgNorm = ( myFlags & EXTRUSION_FLAG_BY_AVG_NORMAL );
+
+  gp_XYZ p = SMESH_TNodeXYZ( srcNode );
+
+  // get normals to faces sharing srcNode
+  vector< gp_XYZ > norms, baryCenters;
+  gp_XYZ norm, avgNorm( 0,0,0 );
+  SMDS_ElemIteratorPtr faceIt = srcNode->GetInverseElementIterator( SMDSAbs_Face );
+  while ( faceIt->more() )
+  {
+    const SMDS_MeshElement* face = faceIt->next();
+    if ( myElemsToUse && !myElemsToUse->count( face ))
+      continue;
+    if ( SMESH_MeshAlgos::FaceNormal( face, norm, /*normalized=*/true ))
+    {
+      norms.push_back( norm );
+      avgNorm += norm;
+      if ( !alongAvgNorm )
+      {
+        gp_XYZ bc(0,0,0);
+        int nbN = 0;
+        for ( SMDS_ElemIteratorPtr nIt = face->nodesIterator(); nIt->more(); ++nbN )
+          bc += SMESH_TNodeXYZ( nIt->next() );
+        baryCenters.push_back( bc / nbN );
+      }
     }
   }
 
-  // create new node and return it
-  const SMDS_MeshNode* NewNode = aMesh->AddNode(x,y,z);
-  //myLastCreatedNodes.Append(NewNode);
-  return NewNode;
+  if ( norms.empty() ) return 0;
+
+  double normSize = avgNorm.Modulus();
+  if ( normSize < std::numeric_limits<double>::min() )
+    return 0;
+
+  if ( myFlags & EXTRUSION_FLAG_BY_AVG_NORMAL ) // extrude along avgNorm
+  {
+    myDir = avgNorm;
+    return makeNodesByDir( mesh, srcNode, newNodes, makeMediumNodes );
+  }
+
+  avgNorm /= normSize;
+
+  int nbNodes = 0;
+  for ( beginStepIter( makeMediumNodes ); moreSteps(); ++nbNodes ) // loop on steps
+  {
+    gp_XYZ pNew = p;
+    double stepSize = nextStep();
+
+    if ( norms.size() > 1 )
+    {
+      for ( size_t iF = 0; iF < norms.size(); ++iF ) // loop on faces
+      {
+        // translate plane of a face
+        baryCenters[ iF ] += norms[ iF ] * stepSize;
+
+        // find point of intersection of the face plane located at baryCenters[ iF ]
+        // and avgNorm located at pNew
+        double d    = -( norms[ iF ] * baryCenters[ iF ]); // d of plane equation ax+by+cz+d=0
+        double dot  = ( norms[ iF ] * avgNorm );
+        if ( dot < std::numeric_limits<double>::min() )
+          dot = stepSize * 1e-3;
+        double step = -( norms[ iF ] * pNew + d ) / dot;
+        pNew += step * avgNorm;
+      }
+    }
+    else
+    {
+      pNew += stepSize * avgNorm;
+    }
+    p = pNew;
+
+    const SMDS_MeshNode * newNode = mesh->AddNode( p.X(), p.Y(), p.Z() );
+    newNodes.push_back( newNode );
+  }
+  return nbNodes;
 }
 
+//=======================================================================
+//function : ExtrusParam::makeNodesByNormal1D
+//purpose  : create nodes for extrusion using normals of edges
+//=======================================================================
+
+int SMESH_MeshEditor::ExtrusParam::
+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 : ExtrusionSweep
@@ -5187,24 +5596,15 @@ const SMDS_MeshNode* SMESH_MeshEditor::CreateNode(const double x,
 //=======================================================================
 
 SMESH_MeshEditor::PGroupIDs
-SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet &   theElems,
+SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet     theElems[2],
                                   const gp_Vec&        theStep,
                                   const int            theNbSteps,
                                   TTElemOfElemListMap& newElemsMap,
-                                  const bool           theMakeGroups,
                                   const int            theFlags,
                                   const double         theTolerance)
 {
-  ExtrusParam aParams;
-  aParams.myDir = gp_Dir(theStep);
-  aParams.myNodes.Clear();
-  aParams.mySteps = new TColStd_HSequenceOfReal;
-  int i;
-  for(i=1; i<=theNbSteps; i++)
-    aParams.mySteps->Append(theStep.Magnitude());
-
-  return
-    ExtrusionSweep(theElems,aParams,newElemsMap,theMakeGroups,theFlags,theTolerance);
+  ExtrusParam aParams( theStep, theNbSteps, theFlags, theTolerance );
+  return ExtrusionSweep( theElems, aParams, newElemsMap );
 }
 
 
@@ -5214,12 +5614,9 @@ SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet &   theElems,
 //=======================================================================
 
 SMESH_MeshEditor::PGroupIDs
-SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet &   theElems,
+SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet     theElemSets[2],
                                   ExtrusParam&         theParams,
-                                  TTElemOfElemListMap& newElemsMap,
-                                  const bool           theMakeGroups,
-                                  const int            theFlags,
-                                  const double         theTolerance)
+                                  TTElemOfElemListMap& newElemsMap)
 {
   myLastCreatedElems.Clear();
   myLastCreatedNodes.Clear();
@@ -5229,7 +5626,9 @@ SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet &   theElems,
 
   SMESHDS_Mesh* aMesh = GetMeshDS();
 
-  int nbsteps = theParams.mySteps->Length();
+  setElemsFirst( theElemSets );
+  const int nbSteps = theParams.NbSteps();
+  theParams.SetElementsToUse( theElemSets[0] );
 
   TNodeOfNodeListMap mapNewNodes;
   //TNodeOfNodeVecMap mapNewNodes;
@@ -5241,91 +5640,75 @@ SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet &   theElems,
                                      myMesh->NbVolumes(ORDER_QUADRATIC) );
   // loop on theElems
   TIDSortedElemSet::iterator itElem;
-  for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
-    // check element type
-    const SMDS_MeshElement* elem = *itElem;
-    if ( !elem  || elem->GetType() == SMDSAbs_Volume )
-      continue;
+  for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
+  {
+    TIDSortedElemSet& theElems = theElemSets[ is2ndSet ];
+    for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
+    {
+      // check element type
+      const SMDS_MeshElement* elem = *itElem;
+      if ( !elem  || elem->GetType() == SMDSAbs_Volume )
+        continue;
 
-    vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
-    newNodesItVec.reserve( elem->NbNodes() );
+      const size_t nbNodes = elem->NbNodes();
+      vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
+      newNodesItVec.reserve( nbNodes );
 
-    // loop on elem nodes
-    SMDS_ElemIteratorPtr itN = elem->nodesIterator();
-    while ( itN->more() )
-    {
-      // check if a node has been already sweeped
-      const SMDS_MeshNode* node = cast2Node( itN->next() );
-      TNodeOfNodeListMap::iterator nIt =
-        mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
-      list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
-      if ( listNewNodes.empty() )
+      // loop on elem nodes
+      SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+      while ( itN->more() )
       {
-        // make new nodes
-
-        // check if we are to create medium nodes between corner ones
-        bool needMediumNodes = false;
-        if ( isQuadraticMesh )
+        // check if a node has been already sweeped
+        const SMDS_MeshNode* node = cast2Node( itN->next() );
+        TNodeOfNodeListMap::iterator nIt =
+          mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
+        list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
+        if ( listNewNodes.empty() )
         {
-          SMDS_ElemIteratorPtr it = node->GetInverseElementIterator();
-          while (it->more() && !needMediumNodes )
-          {
-            const SMDS_MeshElement* invElem = it->next();
-            if ( invElem != elem && !theElems.count( invElem )) continue;
-            needMediumNodes = ( invElem->IsQuadratic() && !invElem->IsMediumNode(node) );
-            if ( !needMediumNodes && invElem->GetEntityType() == SMDSEntity_BiQuad_Quadrangle )
-              needMediumNodes = true;
-          }
-        }
+          // make new nodes
 
-        double coord[] = { node->X(), node->Y(), node->Z() };
-        for ( int i = 0; i < nbsteps; i++ )
-        {
-          if ( needMediumNodes ) // create a medium node
+          // check if we are to create medium nodes between corner ones
+          bool needMediumNodes = false;
+          if ( isQuadraticMesh )
           {
-            double x = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1)/2.;
-            double y = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1)/2.;
-            double z = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1)/2.;
-            if( theFlags & EXTRUSION_FLAG_SEW ) {
-              const SMDS_MeshNode * newNode = CreateNode(x, y, z,
-                                                         theTolerance, theParams.myNodes);
-              listNewNodes.push_back( newNode );
+            SMDS_ElemIteratorPtr it = node->GetInverseElementIterator();
+            while (it->more() && !needMediumNodes )
+            {
+              const SMDS_MeshElement* invElem = it->next();
+              if ( invElem != elem && !theElems.count( invElem )) continue;
+              needMediumNodes = ( invElem->IsQuadratic() && !invElem->IsMediumNode(node) );
+              if ( !needMediumNodes && invElem->GetEntityType() == SMDSEntity_BiQuad_Quadrangle )
+                needMediumNodes = true;
             }
-            else {
-              const SMDS_MeshNode * newNode = aMesh->AddNode(x, y, z);
-              myLastCreatedNodes.Append(newNode);
+          }
+          // create nodes for all steps
+          if ( theParams.MakeNodes( GetMeshDS(), node, listNewNodes, needMediumNodes ))
+          {
+            list<const SMDS_MeshNode*>::iterator newNodesIt = listNewNodes.begin();
+            for ( ; newNodesIt != listNewNodes.end(); ++newNodesIt )
+            {
+              myLastCreatedNodes.Append( *newNodesIt );
               srcNodes.Append( node );
-              listNewNodes.push_back( newNode );
             }
           }
-          // create a corner node
-          coord[0] = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
-          coord[1] = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
-          coord[2] = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
-          if( theFlags & EXTRUSION_FLAG_SEW ) {
-            const SMDS_MeshNode * newNode = CreateNode(coord[0], coord[1], coord[2],
-                                                       theTolerance, theParams.myNodes);
-            listNewNodes.push_back( newNode );
-          }
-          else {
-            const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
-            myLastCreatedNodes.Append(newNode);
-            srcNodes.Append( node );
-            listNewNodes.push_back( newNode );
+          else
+          {
+            break; // newNodesItVec will be shorter than nbNodes
           }
         }
+        newNodesItVec.push_back( nIt );
       }
-      newNodesItVec.push_back( nIt );
+      // make new elements
+      if ( newNodesItVec.size() == nbNodes )
+        sweepElement( elem, newNodesItVec, newElemsMap[elem], nbSteps, srcElems );
     }
-    // make new elements
-    sweepElement( elem, newNodesItVec, newElemsMap[elem], nbsteps, srcElems );
   }
 
-  if( theFlags & EXTRUSION_FLAG_BOUNDARY ) {
-    makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElems, nbsteps, srcElems );
+  if ( theParams.ToMakeBoundary() ) {
+    makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElemSets[0], nbSteps, srcElems );
   }
   PGroupIDs newGroupIDs;
-  if ( theMakeGroups )
+  if ( theParams.ToMakeGroups() )
     newGroupIDs = generateGroups( srcNodes, srcElems, "extruded");
 
   return newGroupIDs;
@@ -5336,7 +5719,7 @@ SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet &   theElems,
 //purpose  :
 //=======================================================================
 SMESH_MeshEditor::Extrusion_Error
-SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
+SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet     theElements[2],
                                        SMESH_subMesh*       theTrack,
                                        const SMDS_MeshNode* theN1,
                                        const bool           theHasAngles,
@@ -5365,7 +5748,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
   TNodeOfNodeListMap mapNewNodes;
 
   // 1. Check data
-  aNbE = theElements.size();
+  aNbE = theElements[0].size() + theElements[1].size();
   // nothing to do
   if ( !aNbE )
     return EXTR_NO_ELEMENTS;
@@ -5508,7 +5891,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
 //purpose  :
 //=======================================================================
 SMESH_MeshEditor::Extrusion_Error
-SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
+SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet     theElements[2],
                                        SMESH_Mesh*          theTrack,
                                        const SMDS_MeshNode* theN1,
                                        const bool           theHasAngles,
@@ -5536,7 +5919,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
   TNodeOfNodeListMap mapNewNodes;
 
   // 1. Check data
-  aNbE = theElements.size();
+  aNbE = theElements[0].size() + theElements[1].size();
   // nothing to do
   if ( !aNbE )
     return EXTR_NO_ELEMENTS;
@@ -5577,7 +5960,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
     }
 
     conn = nbEdgeConnectivity(theN1);
-    if(conn > 2)
+    if( conn != 1 )
       return EXTR_PATH_NOT_EDGE;
 
     aItE = theN1->GetInverseElementIterator();
@@ -5857,7 +6240,7 @@ SMESH_MeshEditor::MakeEdgePathPoints(std::list<double>&                aPrms,
 //purpose  : auxilary for ExtrusionAlongTrack
 //=======================================================================
 SMESH_MeshEditor::Extrusion_Error
-SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet&                 theElements,
+SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet                  theElemSets[2],
                                    list<SMESH_MeshEditor_PathPoint>& fullList,
                                    const bool                        theHasAngles,
                                    list<double>&                     theAngles,
@@ -5895,15 +6278,19 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet&                 theElements
     gp_XYZ aGC( 0.,0.,0. );
     TIDSortedElemSet newNodes;
 
-    itElem = theElements.begin();
-    for ( ; itElem != theElements.end(); itElem++ ) {
-      const SMDS_MeshElement* elem = *itElem;
+    for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
+    {
+      TIDSortedElemSet& theElements = theElemSets[ is2ndSet ];
+      itElem = theElements.begin();
+      for ( ; itElem != theElements.end(); itElem++ ) {
+        const SMDS_MeshElement* elem = *itElem;
 
-      SMDS_ElemIteratorPtr itN = elem->nodesIterator();
-      while ( itN->more() ) {
-        const SMDS_MeshElement* node = itN->next();
-        if ( newNodes.insert( node ).second )
-          aGC += SMESH_TNodeXYZ( node );
+        SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+        while ( itN->more() ) {
+          const SMDS_MeshElement* node = itN->next();
+          if ( newNodes.insert( node ).second )
+            aGC += SMESH_TNodeXYZ( node );
+        }
       }
     }
     aGC /= newNodes.size();
@@ -5913,144 +6300,149 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet&                 theElements
   // 4. Processing the elements
   SMESHDS_Mesh* aMesh = GetMeshDS();
 
-  for ( itElem = theElements.begin(); itElem != theElements.end(); itElem++ ) {
-    // check element type
-    const SMDS_MeshElement* elem = *itElem;
-    SMDSAbs_ElementType   aTypeE = elem->GetType();
-    if ( !elem || ( aTypeE != SMDSAbs_Face && aTypeE != SMDSAbs_Edge ) )
-      continue;
-
-    vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
-    newNodesItVec.reserve( elem->NbNodes() );
+  setElemsFirst( theElemSets );
+  for ( int is2ndSet = 0; is2ndSet < 2; ++is2ndSet )
+  {
+    TIDSortedElemSet& theElements = theElemSets[ is2ndSet ];
+    for ( itElem = theElements.begin(); itElem != theElements.end(); itElem++ ) {
+      // check element type
+      const SMDS_MeshElement* elem = *itElem;
+      SMDSAbs_ElementType   aTypeE = elem->GetType();
+      if ( !elem /*|| ( aTypeE != SMDSAbs_Face && aTypeE != SMDSAbs_Edge )*/ )
+        continue;
 
-    // loop on elem nodes
-    int nodeIndex = -1;
-    SMDS_ElemIteratorPtr itN = elem->nodesIterator();
-    while ( itN->more() )
-    {
-      ++nodeIndex;
-      // check if a node has been already processed
-      const SMDS_MeshNode* node =
-        static_cast<const SMDS_MeshNode*>( itN->next() );
-      TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
-      if ( nIt == mapNewNodes.end() ) {
-        nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
-        list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
+      vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
+      newNodesItVec.reserve( elem->NbNodes() );
 
-        // 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_TNodeXYZ( node );
-
-        const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0];
-        aP0x = aPP0.Pnt();
-        aDT0x= aPP0.Tangent();
-        //cout<<"j = 0   PP: Pnt("<<aP0x.X()<<","<<aP0x.Y()<<","<<aP0x.Z()<<")"<<endl;
-
-        for ( 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 );
-
-          // traslated point
-          aV1x = aV0x.Transformed( aTrsf );
-          aPN1 = aPN0.Transformed( aTrsf );
-
-          // rotation 1 [ T1,T0 ]
-          aAngleT1T0=-aDT1x.Angle( aDT0x );
-          if (fabs(aAngleT1T0) > aTolAng) {
-            aDT1T0=aDT1x^aDT0x;
-            anAxT1T0.SetLocation( aV1x );
-            anAxT1T0.SetDirection( aDT1T0 );
-            aTrsfRotT1T0.SetRotation( anAxT1T0, aAngleT1T0 );
-
-            aPN1 = aPN1.Transformed( aTrsfRotT1T0 );
-          }
+      // loop on elem nodes
+      int nodeIndex = -1;
+      SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+      while ( itN->more() )
+      {
+        ++nodeIndex;
+        // check if a node has been already processed
+        const SMDS_MeshNode* node =
+          static_cast<const SMDS_MeshNode*>( itN->next() );
+        TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
+        if ( nIt == mapNewNodes.end() ) {
+          nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
+          list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
+
+          // make new nodes
+          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_TNodeXYZ( node );
+
+          const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0];
+          aP0x = aPP0.Pnt();
+          aDT0x= aPP0.Tangent();
+          //cout<<"j = 0   PP: Pnt("<<aP0x.X()<<","<<aP0x.Y()<<","<<aP0x.Z()<<")"<<endl;
+
+          for ( 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 );
+
+            // traslated point
+            aV1x = aV0x.Transformed( aTrsf );
+            aPN1 = aPN0.Transformed( aTrsf );
+
+            // rotation 1 [ T1,T0 ]
+            aAngleT1T0=-aDT1x.Angle( aDT0x );
+            if (fabs(aAngleT1T0) > aTolAng) {
+              aDT1T0=aDT1x^aDT0x;
+              anAxT1T0.SetLocation( aV1x );
+              anAxT1T0.SetDirection( aDT1T0 );
+              aTrsfRotT1T0.SetRotation( anAxT1T0, aAngleT1T0 );
+
+              aPN1 = aPN1.Transformed( aTrsfRotT1T0 );
+            }
 
-          // rotation 2
-          if ( theHasAngles ) {
-            anAx1.SetLocation( aV1x );
-            anAx1.SetDirection( aDT1x );
-            aTrsfRot.SetRotation( anAx1, aAngle1x );
+            // rotation 2
+            if ( theHasAngles ) {
+              anAx1.SetLocation( aV1x );
+              anAx1.SetDirection( aDT1x );
+              aTrsfRot.SetRotation( anAx1, aAngle1x );
 
-            aPN1 = aPN1.Transformed( aTrsfRot );
-          }
+              aPN1 = aPN1.Transformed( aTrsfRot );
+            }
 
-          // make new node
-          //MESSAGE("elem->IsQuadratic " << elem->IsQuadratic() << " " << elem->IsMediumNode(node));
-          if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
-            // create additional node
-            double x = ( aPN1.X() + aPN0.X() )/2.;
-            double y = ( aPN1.Y() + aPN0.Y() )/2.;
-            double z = ( aPN1.Z() + aPN0.Z() )/2.;
-            const SMDS_MeshNode* newNode = aMesh->AddNode(x,y,z);
+            // make new node
+            //MESSAGE("elem->IsQuadratic " << elem->IsQuadratic() << " " << elem->IsMediumNode(node));
+            if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
+              // create additional node
+              double x = ( aPN1.X() + aPN0.X() )/2.;
+              double y = ( aPN1.Y() + aPN0.Y() )/2.;
+              double z = ( aPN1.Z() + aPN0.Z() )/2.;
+              const SMDS_MeshNode* newNode = aMesh->AddNode(x,y,z);
+              myLastCreatedNodes.Append(newNode);
+              srcNodes.Append( node );
+              listNewNodes.push_back( newNode );
+            }
+            const SMDS_MeshNode* newNode = aMesh->AddNode( aPN1.X(), aPN1.Y(), aPN1.Z() );
             myLastCreatedNodes.Append(newNode);
             srcNodes.Append( node );
             listNewNodes.push_back( newNode );
-          }
-          const SMDS_MeshNode* newNode = aMesh->AddNode( aPN1.X(), aPN1.Y(), aPN1.Z() );
-          myLastCreatedNodes.Append(newNode);
-          srcNodes.Append( node );
-          listNewNodes.push_back( newNode );
 
-          aPN0 = aPN1;
-          aP0x = aP1x;
-          aV0x = aV1x;
-          aDT0x = aDT1x;
+            aPN0 = aPN1;
+            aP0x = aP1x;
+            aV0x = aV1x;
+            aDT0x = aDT1x;
+          }
         }
-      }
 
-      else {
-        // if current elem is quadratic and current node is not medium
-        // we have to check - may be it is needed to insert additional nodes
-        if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
-          list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
-          if(listNewNodes.size()==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.Append( node );
-              myLastCreatedNodes.Append(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]);
+        else {
+          // if current elem is quadratic and current node is not medium
+          // we have to check - may be it is needed to insert additional nodes
+          if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
+            list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
+            if(listNewNodes.size()==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.Append( node );
+                myLastCreatedNodes.Append(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]);
+              }
             }
           }
         }
-      }
 
-      newNodesItVec.push_back( nIt );
+        newNodesItVec.push_back( nIt );
+      }
+      // make new elements
+      //sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem],
+      //              newNodesItVec[0]->second.size(), myLastCreatedElems );
+      sweepElement( elem, newNodesItVec, newElemsMap[elem], aNbTP-1, srcElems );
     }
-    // make new elements
-    //sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem],
-    //              newNodesItVec[0]->second.size(), myLastCreatedElems );
-    sweepElement( elem, newNodesItVec, newElemsMap[elem], aNbTP-1, srcElems );
   }
 
-  makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElements, aNbTP-1, srcElems );
+  makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElemSets[0], aNbTP-1, srcElems );
 
   if ( theMakeGroups )
     generateGroups( srcNodes, srcElems, "extruded");
@@ -9988,7 +10380,7 @@ SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
 
 void SMESH_MeshEditor::DoubleElements( const TIDSortedElemSet& theElements )
 {
-  CrearLastCreated();
+  ClearLastCreated();
   SMESHDS_Mesh* mesh = GetMeshDS();
 
   // get an element type and an iterator over elements
@@ -10147,6 +10539,7 @@ bool SMESH_MeshEditor::doubleNodes( SMESHDS_Mesh*     theMeshDS,
       {
         // duplicate node
         aNewNode = theMeshDS->AddNode( aCurrNode->X(), aCurrNode->Y(), aCurrNode->Z() );
+        copyPosition( aCurrNode, aNewNode );
         theNodeNodeMap[ aCurrNode ] = aNewNode;
         myLastCreatedNodes.Append( aNewNode );
       }
@@ -10159,10 +10552,8 @@ bool SMESH_MeshEditor::doubleNodes( SMESHDS_Mesh*     theMeshDS,
     if ( theIsDoubleElem )
       AddElement(newNodes, anElem->GetType(), anElem->IsPoly());
     else
-      {
-      MESSAGE("ChangeElementNodes");
       theMeshDS->ChangeElementNodes( anElem, &newNodes[ 0 ], anElem->NbNodes() );
-      }
+
     res = true;
   }
   return res;
@@ -10173,8 +10564,8 @@ bool SMESH_MeshEditor::doubleNodes( SMESHDS_Mesh*     theMeshDS,
   \brief Creates a hole in a mesh by doubling the nodes of some particular elements
   \param theNodes - identifiers of nodes to be doubled
   \param theModifiedElems - identifiers of elements to be updated by the new (doubled)
-         nodes. If list of element identifiers is empty then nodes are doubled but
-         they not assigned to elements
+  nodes. If list of element identifiers is empty then nodes are doubled but
+  they not assigned to elements
   \return TRUE if operation has been completed successfully, FALSE otherwise
 */
 //================================================================================
@@ -10210,6 +10601,7 @@ bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes,
     const SMDS_MeshNode* aNewNode = aMeshDS->AddNode( aNode->X(), aNode->Y(), aNode->Z() );
     if ( aNewNode )
     {
+      copyPosition( aNode, aNewNode );
       anOldNodeToNewNode[ aNode ] = aNewNode;
       myLastCreatedNodes.Append( aNewNode );
     }
@@ -10311,15 +10703,12 @@ namespace {
     }
     void Perform(const gp_Pnt& aPnt, double theTol)
     {
+      theTol *= theTol;
       _state = TopAbs_OUT;
       _extremum.Perform(aPnt);
       if ( _extremum.IsDone() )
         for ( int iSol = 1; iSol <= _extremum.NbExt() && _state == TopAbs_OUT; ++iSol)
-#if OCC_VERSION_LARGE > 0x06040000 // Porting to OCCT6.5.1
           _state = ( _extremum.SquareDistance(iSol) <= theTol ? TopAbs_IN : TopAbs_OUT );
-#else
-          _state = ( _extremum.Value(iSol) <= theTol ? TopAbs_IN : TopAbs_OUT );
-#endif
     }
     TopAbs_State State() const
     {
@@ -10831,6 +11220,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
                         }
                       double *coords = grid->GetPoint(oldId);
                       SMDS_MeshNode *newNode = meshDS->AddNode(coords[0], coords[1], coords[2]);
+                      copyPosition( meshDS->FindNodeVtk( oldId ), newNode );
                       int newId = newNode->getVtkId();
                       nodeDomains[oldId][idom] = newId; // cloned node for other domains
                       //MESSAGE("-+-+-c     oldNode " << oldId << " domain " << idomain << " newNode " << newId << " domain " << idom << " size=" <<nodeDomains[oldId].size());
@@ -11169,6 +11559,14 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
         }
     }
 
+  // Remove empty groups (issue 0022812)
+  std::map<std::string, SMESH_Group*>::iterator name_group = mapOfJunctionGroups.begin();
+  for ( ; name_group != mapOfJunctionGroups.end(); ++name_group )
+  {
+    if ( name_group->second && name_group->second->GetGroupDS()->IsEmpty() )
+      myMesh->RemoveGroup( name_group->second->GetGroupDS()->GetID() );
+  }
+
   meshDS->CleanDownWardConnectivity(); // Mesh has been modified, downward connectivity is no more usable, free memory
   grid->BuildLinks();
 
@@ -11236,6 +11634,7 @@ bool SMESH_MeshEditor::CreateFlatElementsOnFacesGroups(const std::vector<TIDSort
               if (!clonedNodes.count(node))
                 {
                   clone = meshDS->AddNode(node->X(), node->Y(), node->Z());
+                  copyPosition( node, clone );
                   clonedNodes[node] = clone;
                 }
               else
@@ -11252,6 +11651,7 @@ bool SMESH_MeshEditor::CreateFlatElementsOnFacesGroups(const std::vector<TIDSort
                   if (!intermediateNodes.count(node))
                     {
                       inter = meshDS->AddNode(node->X(), node->Y(), node->Z());
+                      copyPosition( node, inter );
                       intermediateNodes[node] = inter;
                     }
                   else
@@ -12196,3 +12596,45 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
   }
   return nbAddedBnd;
 }
+
+//================================================================================
+/*!
+ * \brief Copy node position and set \a to node on the same geometry
+ */
+//================================================================================
+
+void SMESH_MeshEditor::copyPosition( const SMDS_MeshNode* from,
+                                     const SMDS_MeshNode* to )
+{
+  if ( !from || !to ) return;
+
+  SMDS_PositionPtr pos = from->GetPosition();
+  if ( !pos || from->getshapeId() < 1 ) return;
+
+  switch ( pos->GetTypeOfPosition() )
+  {
+  case SMDS_TOP_3DSPACE: break;
+
+  case SMDS_TOP_FACE:
+  {
+    const SMDS_FacePosition* fPos = static_cast< const SMDS_FacePosition* >( pos );
+    GetMeshDS()->SetNodeOnFace( to, from->getshapeId(),
+                                fPos->GetUParameter(), fPos->GetVParameter() );
+    break;
+  }
+  case SMDS_TOP_EDGE:
+  {
+    // WARNING: it is dangerous to set equal nodes on one EDGE!!!!!!!!
+    const SMDS_EdgePosition* ePos = static_cast< const SMDS_EdgePosition* >( pos );
+    GetMeshDS()->SetNodeOnEdge( to, from->getshapeId(), ePos->GetUParameter() );
+    break;
+  }
+  case SMDS_TOP_VERTEX:
+  {
+    GetMeshDS()->SetNodeOnVertex( to, from->getshapeId() );
+    break;
+  }
+  case SMDS_TOP_UNSPEC:
+  default:;
+  }
+}