Salome HOME
Regression of smesh/mesh_Projection_2D_01/B7
[modules/smesh.git] / src / SMESH / SMESH_MeshEditor.cxx
index 842b177507e10c3f190e4a1d626a26ca6c4a1bff..b45e0117bb235707b4499f60e482f6619c487794 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2013  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
@@ -6,7 +6,7 @@
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 // License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
 //
 // This library is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
@@ -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();
@@ -518,14 +519,12 @@ int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem)
   }
   else
   {
-    const map<int,SMESHDS_SubMesh*>& id2sm = GetMeshDS()->SubMeshes();
-    map<int,SMESHDS_SubMesh*>::const_iterator id_sm = id2sm.begin();
-    for ( ; id_sm != id2sm.end(); ++id_sm )
-      if ( id_sm->second->Contains( theElem ))
-        return id_sm->first;
+    SMESHDS_SubMeshIteratorPtr smIt = GetMeshDS()->SubMeshes();
+    while ( const SMESHDS_SubMesh* sm = smIt->next() )
+      if ( sm->Contains( theElem ))
+        return sm->GetID();
   }
 
-  //MESSAGE ("::FindShape() - SHAPE NOT FOUND")
   return 0;
 }
 
@@ -1174,7 +1173,7 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
     avoidSet.clear();
     avoidSet.insert(theFace);
 
-    NLink link( theFace->GetNode( 0 ), 0 );
+    NLink link( theFace->GetNode( 0 ), (SMDS_MeshNode *) 0 );
 
     const int nbNodes = theFace->NbCornerNodes();
     for ( int i = 0; i < nbNodes; ++i ) // loop on links of theFace
@@ -1245,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  :
@@ -1603,44 +1688,110 @@ namespace
   const int* thePentaTo3[6] = { thePentaTo3_1, thePentaTo3_2, thePentaTo3_3,
                                 thePentaTo3_4, thePentaTo3_5, thePentaTo3_6 };
 
+  // Methods of splitting hexahedron into prisms
+
+  const int theHexTo4Prisms_BT[6*4+1] = // bottom-top
+    {
+      0, 1, 8, 4, 5, 9,    1, 2, 8, 5, 6, 9,    2, 3, 8, 6, 7, 9,   3, 0, 8, 7, 4, 9,    -1
+    };
+  const int theHexTo4Prisms_LR[6*4+1] = // left-right
+    {
+      1, 0, 8, 2, 3, 9,    0, 4, 8, 3, 7, 9,    4, 5, 8, 7, 6, 9,   5, 1, 8, 6, 2, 9,    -1
+    };
+  const int theHexTo4Prisms_FB[6*4+1] = // front-back
+    {
+      0, 3, 9, 1, 2, 8,    3, 7, 9, 2, 6, 8,    7, 4, 9, 6, 5, 8,   4, 0, 9, 5, 1, 8,    -1
+    };
+
+  const int theHexTo2Prisms_BT_1[6*2+1] =
+    {
+      0, 1, 3, 4, 5, 7,    1, 2, 3, 5, 6, 7,   -1
+    };
+  const int theHexTo2Prisms_BT_2[6*2+1] =
+    {
+      0, 1, 2, 4, 5, 6,    0, 2, 3, 4, 6, 7,   -1
+    };
+  const int* theHexTo2Prisms_BT[2] = { theHexTo2Prisms_BT_1, theHexTo2Prisms_BT_2 };
+
+  const int theHexTo2Prisms_LR_1[6*2+1] =
+    {
+      1, 0, 4, 2, 3, 7,    1, 4, 5, 2, 7, 6,   -1
+    };
+  const int theHexTo2Prisms_LR_2[6*2+1] =
+    {
+      1, 0, 4, 2, 3, 7,    1, 4, 5, 2, 7, 6,   -1
+    };
+  const int* theHexTo2Prisms_LR[2] = { theHexTo2Prisms_LR_1, theHexTo2Prisms_LR_2 };
+
+  const int theHexTo2Prisms_FB_1[6*2+1] =
+    {
+      0, 3, 4, 1, 2, 5,    3, 7, 4, 2, 6, 5,   -1
+    };
+  const int theHexTo2Prisms_FB_2[6*2+1] =
+    {
+      0, 3, 7, 1, 2, 7,    0, 7, 4, 1, 6, 5,   -1
+    };
+  const int* theHexTo2Prisms_FB[2] = { theHexTo2Prisms_FB_1, theHexTo2Prisms_FB_2 };
+
+
   struct TTriangleFacet //!< stores indices of three nodes of tetra facet
   {
     int _n1, _n2, _n3;
     TTriangleFacet(int n1, int n2, int n3): _n1(n1), _n2(n2), _n3(n3) {}
     bool contains(int n) const { return ( n == _n1 || n == _n2 || n == _n3 ); }
-    bool hasAdjacentTetra( const SMDS_MeshElement* elem ) const;
+    bool hasAdjacentVol( const SMDS_MeshElement*    elem,
+                         const SMDSAbs_GeometryType geom = SMDSGeom_TETRA) const;
   };
   struct TSplitMethod
   {
-    int        _nbTetra;
+    int        _nbSplits;
+    int        _nbCorners;
     const int* _connectivity; //!< foursomes of tetra connectivy finished by -1
     bool       _baryNode;     //!< additional node is to be created at cell barycenter
     bool       _ownConn;      //!< to delete _connectivity in destructor
     map<int, const SMDS_MeshNode*> _faceBaryNode; //!< map face index to node at BC of face
 
     TSplitMethod( int nbTet=0, const int* conn=0, bool addNode=false)
-      : _nbTetra(nbTet), _connectivity(conn), _baryNode(addNode), _ownConn(false) {}
+      : _nbSplits(nbTet), _nbCorners(4), _connectivity(conn), _baryNode(addNode), _ownConn(false) {}
     ~TSplitMethod() { if ( _ownConn ) delete [] _connectivity; _connectivity = 0; }
     bool hasFacet( const TTriangleFacet& facet ) const
     {
-      const int* tetConn = _connectivity;
-      for ( ; tetConn[0] >= 0; tetConn += 4 )
-        if (( facet.contains( tetConn[0] ) +
-              facet.contains( tetConn[1] ) +
-              facet.contains( tetConn[2] ) +
-              facet.contains( tetConn[3] )) == 3 )
-          return true;
+      if ( _nbCorners == 4 )
+      {
+        const int* tetConn = _connectivity;
+        for ( ; tetConn[0] >= 0; tetConn += 4 )
+          if (( facet.contains( tetConn[0] ) +
+                facet.contains( tetConn[1] ) +
+                facet.contains( tetConn[2] ) +
+                facet.contains( tetConn[3] )) == 3 )
+            return true;
+      }
+      else // prism, _nbCorners == 6
+      {
+        const int* prismConn = _connectivity;
+        for ( ; prismConn[0] >= 0; prismConn += 6 )
+        {
+          if (( facet.contains( prismConn[0] ) &&
+                facet.contains( prismConn[1] ) &&
+                facet.contains( prismConn[2] ))
+              ||
+              ( facet.contains( prismConn[3] ) &&
+                facet.contains( prismConn[4] ) &&
+                facet.contains( prismConn[5] )))
+            return true;
+        }
+      }
       return false;
     }
   };
 
   //=======================================================================
   /*!
-   * \brief return TSplitMethod for the given element
+   * \brief return TSplitMethod for the given element to split into tetrahedra
    */
   //=======================================================================
 
-  TSplitMethod getSplitMethod( SMDS_VolumeTool& vol, const int theMethodFlags)
+  TSplitMethod getTetraSplitMethod( SMDS_VolumeTool& vol, const int theMethodFlags)
   {
     const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
 
@@ -1665,8 +1816,8 @@ namespace
       {
         TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] );
         TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] );
-        if      ( t012.hasAdjacentTetra( vol.Element() )) triaSplits.push_back( t012 );
-        else if ( t123.hasAdjacentTetra( vol.Element() )) triaSplits.push_back( t123 );
+        if      ( t012.hasAdjacentVol( vol.Element() )) triaSplits.push_back( t012 );
+        else if ( t123.hasAdjacentVol( vol.Element() )) triaSplits.push_back( t123 );
       }
       else
       {
@@ -1679,7 +1830,7 @@ namespace
           TTriangleFacet t023( nInd[ iQ * ( iCom             )],
                                nInd[ iQ * ( (iCom+2)%nbNodes )],
                                nInd[ iQ * ( (iCom+3)%nbNodes )]);
-          if ( t012.hasAdjacentTetra( vol.Element() ) && t023.hasAdjacentTetra( vol.Element() ))
+          if ( t012.hasAdjacentVol( vol.Element() ) && t023.hasAdjacentVol( vol.Element() ))
           {
             triaSplits.push_back( t012 );
             triaSplits.push_back( t023 );
@@ -1719,12 +1870,12 @@ namespace
       default:
         nbVariants = 0;
       }
-      for ( int variant = 0; variant < nbVariants && method._nbTetra == 0; ++variant )
+      for ( int variant = 0; variant < nbVariants && method._nbSplits == 0; ++variant )
       {
         // check method compliancy with adjacent tetras,
         // all found splits must be among facets of tetras described by this method
         method = TSplitMethod( nbTet, connVariants[variant] );
-        if ( hasAdjacentSplits && method._nbTetra > 0 )
+        if ( hasAdjacentSplits && method._nbSplits > 0 )
         {
           bool facetCreated = true;
           for ( int iF = 0; facetCreated && iF < triaSplitsByFace.size(); ++iF )
@@ -1738,7 +1889,7 @@ namespace
         }
       }
     }
-    if ( method._nbTetra < 1 )
+    if ( method._nbSplits < 1 )
     {
       // No standard method is applicable, use a generic solution:
       // each facet of a volume is split into triangles and
@@ -1832,7 +1983,7 @@ namespace
             connectivity[ connSize++ ] = baryCenInd;
           }
         }
-        method._nbTetra += nbTet;
+        method._nbSplits += nbTet;
 
       } // loop on volume faces
 
@@ -1842,13 +1993,132 @@ namespace
 
     return method;
   }
+  //=======================================================================
+  /*!
+   * \brief return TSplitMethod to split haxhedron into prisms
+   */
+  //=======================================================================
+
+  TSplitMethod getPrismSplitMethod( SMDS_VolumeTool& vol,
+                                    const int        methodFlags,
+                                    const int        facetToSplit)
+  {
+    // 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]
+
+    if ( methodFlags == SMESH_MeshEditor::HEXA_TO_4_PRISMS )
+    {
+      static TSplitMethod to4methods[4]; // order BT, LR, FB
+      if ( to4methods[iF]._nbSplits == 0 )
+      {
+        switch ( iF ) {
+        case 0:
+          to4methods[iF]._connectivity = theHexTo4Prisms_BT;
+          to4methods[iF]._faceBaryNode[ 0 ] = 0;
+          to4methods[iF]._faceBaryNode[ 1 ] = 0;
+          break;
+        case 1:
+          to4methods[iF]._connectivity = theHexTo4Prisms_LR;
+          to4methods[iF]._faceBaryNode[ 2 ] = 0;
+          to4methods[iF]._faceBaryNode[ 4 ] = 0;
+          break;
+        case 2:
+          to4methods[iF]._connectivity = theHexTo4Prisms_FB;
+          to4methods[iF]._faceBaryNode[ 3 ] = 0;
+          to4methods[iF]._faceBaryNode[ 5 ] = 0;
+          break;
+        default: return to4methods[3];
+        }
+        to4methods[iF]._nbSplits  = 4;
+        to4methods[iF]._nbCorners = 6;
+      }
+      return to4methods[iF];
+    }
+    // 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** connVariants = 0;
+    switch ( iF ) {
+    case 0: connVariants = theHexTo2Prisms_BT; break;
+    case 1: connVariants = theHexTo2Prisms_LR; break;
+    case 2: connVariants = theHexTo2Prisms_FB; break;
+    default: return method;
+    }
+
+    // look for prisms adjacent via facetToSplit and an opposite one
+    for ( int is2nd = 0; is2nd < 2; ++is2nd )
+    {
+      int iFacet = is2nd ? vol.GetOppFaceIndexOfHex( facetToSplit ) : facetToSplit;
+      int nbNodes = vol.NbFaceNodes( iFacet ) / iQ;
+      if ( nbNodes != 4 ) return method;
+
+      const int* nInd = vol.GetFaceNodesIndices( iFacet );
+      TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] );
+      TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] );
+      TTriangleFacet* t;
+      if      ( t012.hasAdjacentVol( vol.Element(), SMDSGeom_PENTA ))
+        t = &t012;
+      else if ( t123.hasAdjacentVol( vol.Element(), SMDSGeom_PENTA ))
+        t = &t123;
+      else
+        continue;
+
+      // there are adjacent prism
+      for ( int variant = 0; variant < nbVariants; ++variant )
+      {
+        // check method compliancy with adjacent prisms,
+        // the found prism facets must be among facets of prisms described by current method
+        method._nbSplits     = nbSplits;
+        method._nbCorners    = 6;
+        method._connectivity = connVariants[ variant ];
+        if ( method.hasFacet( *t ))
+          return method;
+      }
+    }
+
+    // No adjacent prisms. Select a variant with a best aspect ratio.
+
+    double badness[2] = { 0, 0 };
+    static SMESH::Controls::NumericalFunctorPtr aspectRatio( new SMESH::Controls::AspectRatio);
+    const SMDS_MeshNode** nodes = vol.GetNodes();
+    for ( int variant = 0; variant < nbVariants; ++variant )
+      for ( int is2nd = 0; is2nd < 2; ++is2nd )
+      {
+        int iFacet = is2nd ? vol.GetOppFaceIndexOfHex( facetToSplit ) : facetToSplit;
+        const int*             nInd = vol.GetFaceNodesIndices( iFacet );
+
+        method._connectivity = connVariants[ variant ];
+        TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] );
+        TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] );
+        TTriangleFacet* t = ( method.hasFacet( t012 )) ? & t012 : & t123;
+
+        SMDS_FaceOfNodes tria ( nodes[ t->_n1 ],
+                                nodes[ t->_n2 ],
+                                nodes[ t->_n3 ] );
+        badness[ variant ] += getBadRate( &tria, aspectRatio );
+      }
+    const int iBetter = ( badness[1] < badness[0] && badness[0]-badness[1] > 0.1 * badness[0] );
+
+    method._nbSplits     = nbSplits;
+    method._nbCorners    = 6;
+    method._connectivity = connVariants[ iBetter ];
+
+    return method;
+  }
+
   //================================================================================
   /*!
    * \brief Check if there is a tetraherdon adjacent to the given element via this facet
    */
   //================================================================================
 
-  bool TTriangleFacet::hasAdjacentTetra( const SMDS_MeshElement* elem ) const
+  bool TTriangleFacet::hasAdjacentVol( const SMDS_MeshElement*    elem,
+                                       const SMDSAbs_GeometryType geom ) const
   {
     // find the tetrahedron including the three nodes of facet
     const SMDS_MeshNode* n1 = elem->GetNode(_n1);
@@ -1858,16 +2128,16 @@ namespace
     while ( volIt1->more() )
     {
       const SMDS_MeshElement* v = volIt1->next();
-      SMDSAbs_EntityType type = v->GetEntityType();
-      if ( type != SMDSEntity_Tetra && type != SMDSEntity_Quad_Tetra )
+      if ( v->GetGeomType() != geom )
         continue;
-      if ( type == SMDSEntity_Quad_Tetra && v->GetNodeIndex( n1 ) > 3 )
+      const int lastCornerInd = v->NbCornerNodes() - 1;
+      if ( v->IsQuadratic() && v->GetNodeIndex( n1 ) > lastCornerInd )
         continue; // medium node not allowed
       const int ind2 = v->GetNodeIndex( n2 );
-      if ( ind2 < 0 || 3 < ind2 )
+      if ( ind2 < 0 || lastCornerInd < ind2 )
         continue;
       const int ind3 = v->GetNodeIndex( n3 );
-      if ( ind3 < 0 || 3 < ind3 )
+      if ( ind3 < 0 || lastCornerInd < ind3 )
         continue;
       return true;
     }
@@ -1900,19 +2170,19 @@ namespace
 } // namespace
 
 //=======================================================================
-//function : SplitVolumesIntoTetra
-//purpose  : Split volume elements into tetrahedra.
+//function : SplitVolumes
+//purpose  : Split volume elements into tetrahedra or prisms.
+//           If facet ID < 0, element is split into tetrahedra,
+//           else a hexahedron is split into prisms so that the given facet is
+//           split into triangles
 //=======================================================================
 
-void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
-                                              const int                theMethodFlags)
+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());
+  SMESH_MesherHelper helper( *GetMesh()), fHelper(*GetMesh());
+  fHelper.ToFixNodeParameters( true );
 
   SMESHDS_SubMesh* subMesh = 0;//GetMeshDS()->MeshElements(1);
   SMESHDS_SubMesh* fSubMesh = 0;//subMesh;
@@ -1923,29 +2193,33 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
   map< TVolumeFaceKey, const SMDS_MeshNode* > volFace2BaryNode;
   double bc[3];
 
-  TIDSortedElemSet::const_iterator elem = theElems.begin();
-  for ( ; elem != theElems.end(); ++elem )
+  TFacetOfElem::const_iterator elem2facet = theElems.begin();
+  for ( ; elem2facet != theElems.end(); ++elem2facet )
   {
-    if ( (*elem)->GetType() != SMDSAbs_Volume )
+    const SMDS_MeshElement* elem = elem2facet->first;
+    const int       facetToSplit = elem2facet->second;
+    if ( elem->GetType() != SMDSAbs_Volume )
       continue;
-    SMDSAbs_EntityType geomType = (*elem)->GetEntityType();
+    const SMDSAbs_EntityType geomType = elem->GetEntityType();
     if ( geomType == SMDSEntity_Tetra || geomType == SMDSEntity_Quad_Tetra )
       continue;
 
-    if ( !volTool.Set( *elem, /*ignoreCentralNodes=*/false )) continue; // strange...
+    if ( !volTool.Set( elem, /*ignoreCentralNodes=*/false )) continue; // strange...
 
-    TSplitMethod splitMethod = getSplitMethod( volTool, theMethodFlags );
-    if ( splitMethod._nbTetra < 1 ) continue;
+    TSplitMethod splitMethod = ( facetToSplit < 0  ?
+                                 getTetraSplitMethod( volTool, theMethodFlags ) :
+                                 getPrismSplitMethod( volTool, theMethodFlags, facetToSplit ));
+    if ( splitMethod._nbSplits < 1 ) continue;
 
     // find submesh to add new tetras to
-    if ( !subMesh || !subMesh->Contains( *elem ))
+    if ( !subMesh || !subMesh->Contains( elem ))
     {
-      int shapeID = FindShape( *elem );
+      int shapeID = FindShape( elem );
       helper.SetSubShape( shapeID ); // helper will add tetras to the found submesh
       subMesh = GetMeshDS()->MeshElements( shapeID );
     }
     int iQ;
-    if ( (*elem)->IsQuadratic() )
+    if ( elem->IsQuadratic() )
     {
       iQ = 2;
       // add quadratic links to the helper
@@ -1963,7 +2237,8 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
       iQ = 1;
       helper.SetIsQuadratic( false );
     }
-    vector<const SMDS_MeshNode*> nodes( (*elem)->begin_nodes(), (*elem)->end_nodes() );
+    vector<const SMDS_MeshNode*> nodes( volTool.GetNodes(),
+                                        volTool.GetNodes() + elem->NbNodes() );
     helper.SetElementsOnShape( true );
     if ( splitMethod._baryNode )
     {
@@ -1991,16 +2266,25 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
       }
     }
 
-    // make tetras
-    vector<const SMDS_MeshElement* > tetras( splitMethod._nbTetra ); // splits of a volume
-    const int* tetConn = splitMethod._connectivity;
-    for ( int i = 0; i < splitMethod._nbTetra; ++i, tetConn += 4 )
-      newElems.Append( tetras[ i ] = helper.AddVolume( nodes[ tetConn[0] ],
-                                                       nodes[ tetConn[1] ],
-                                                       nodes[ tetConn[2] ],
-                                                       nodes[ tetConn[3] ]));
+    // make new volumes
+    vector<const SMDS_MeshElement* > splitVols( splitMethod._nbSplits ); // splits of a volume
+    const int* volConn = splitMethod._connectivity;
+    if ( splitMethod._nbCorners == 4 ) // tetra
+      for ( int i = 0; i < splitMethod._nbSplits; ++i, volConn += splitMethod._nbCorners )
+        newElems.Append( splitVols[ i ] = helper.AddVolume( nodes[ volConn[0] ],
+                                                            nodes[ volConn[1] ],
+                                                            nodes[ volConn[2] ],
+                                                            nodes[ volConn[3] ]));
+    else // prisms
+      for ( int i = 0; i < splitMethod._nbSplits; ++i, volConn += splitMethod._nbCorners )
+        newElems.Append( splitVols[ i ] = helper.AddVolume( nodes[ volConn[0] ],
+                                                            nodes[ volConn[1] ],
+                                                            nodes[ volConn[2] ],
+                                                            nodes[ volConn[3] ],
+                                                            nodes[ volConn[4] ],
+                                                            nodes[ volConn[5] ]));
 
-    ReplaceElemInGroups( *elem, tetras, GetMeshDS() );
+    ReplaceElemInGroups( elem, splitVols, GetMeshDS() );
 
     // Split faces on sides of the split volume
 
@@ -2029,17 +2313,37 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
         map<int, const SMDS_MeshNode*>::iterator iF_n = splitMethod._faceBaryNode.find(iF);
         if ( iF_n != splitMethod._faceBaryNode.end() )
         {
+          const SMDS_MeshNode *baryNode = iF_n->second;
           for ( int iN = 0; iN < nbNodes*iQ; iN += iQ )
           {
             const SMDS_MeshNode* n1 = fNodes[iN];
             const SMDS_MeshNode *n2 = fNodes[(iN+iQ)%(nbNodes*iQ)];
-            const SMDS_MeshNode *n3 = iF_n->second;
+            const SMDS_MeshNode *n3 = baryNode;
             if ( !volTool.IsFaceExternal( iF ))
               swap( n2, n3 );
             triangles.push_back( helper.AddFace( n1,n2,n3 ));
-
-            if ( fSubMesh && n3->getshapeId() < 1 )
-              fSubMesh->AddNode( n3 );
+          }
+          if ( fSubMesh ) // update position of the bary node on geometry
+          {
+            if ( subMesh )
+              subMesh->RemoveNode( baryNode, false );
+            GetMeshDS()->SetNodeOnFace( baryNode, fSubMesh->GetID() );
+            const TopoDS_Shape& s = GetMeshDS()->IndexToShape( fSubMesh->GetID() );
+            if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
+            {
+              fHelper.SetSubShape( s );
+              gp_XY uv( 1e100, 1e100 );
+              double distXYZ[4];
+              if ( !fHelper.CheckNodeUV( TopoDS::Face( s ), baryNode,
+                                        uv, /*tol=*/1e-7, /*force=*/true, distXYZ ) &&
+                   uv.X() < 1e100 )
+              {
+                // node is too far from the surface
+                GetMeshDS()->MoveNode( baryNode, distXYZ[1], distXYZ[2], distXYZ[3] );
+                const_cast<SMDS_MeshNode*>( baryNode )->SetPosition
+                  ( SMDS_PositionPtr( new SMDS_FacePosition( uv.X(), uv.Y() )));
+              }
+            }
           }
         }
         else
@@ -2069,6 +2373,8 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
             }
           }
           list< TTriangleFacet >::iterator facet = facets.begin();
+          if ( facet == facets.end() )
+            break;
           for ( ; facet != facets.end(); ++facet )
           {
             if ( !volTool.IsFaceExternal( iF ))
@@ -2087,11 +2393,11 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
         }
         ReplaceElemInGroups( face, triangles, GetMeshDS() );
         GetMeshDS()->RemoveFreeElement( face, fSubMesh, /*fromGroups=*/false );
-      }
 
+      } // while a face based on facet nodes exists
     } // loop on volume faces to split them into triangles
 
-    GetMeshDS()->RemoveFreeElement( *elem, subMesh, /*fromGroups=*/false );
+    GetMeshDS()->RemoveFreeElement( elem, subMesh, /*fromGroups=*/false );
 
     if ( geomType == SMDSEntity_TriQuad_Hexa )
     {
@@ -2106,6 +2412,198 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
   myLastCreatedElems = newElems;
 }
 
+//=======================================================================
+//function : GetHexaFacetsToSplit
+//purpose  : For hexahedra that will be split into prisms, finds facets to
+//           split into triangles. Only hexahedra adjacent to the one closest
+//           to theFacetNormal.Location() are returned.
+//param [in,out] theHexas - the hexahedra
+//param [in]     theFacetNormal - facet normal
+//param [out]    theFacets - the hexahedra and found facet IDs
+//=======================================================================
+
+void SMESH_MeshEditor::GetHexaFacetsToSplit( TIDSortedElemSet& theHexas,
+                                             const gp_Ax1&     theFacetNormal,
+                                             TFacetOfElem &    theFacets)
+{
+  #define THIS_METHOD "SMESH_MeshEditor::GetHexaFacetsToSplit(): "
+
+  // Find a hexa closest to the location of theFacetNormal
+
+  const SMDS_MeshElement* startHex;
+  {
+    // get SMDS_ElemIteratorPtr on theHexas
+    typedef const SMDS_MeshElement*                                      TValue;
+    typedef TIDSortedElemSet::iterator                                   TSetIterator;
+    typedef SMDS::SimpleAccessor<TValue,TSetIterator>                    TAccesor;
+    typedef SMDS_MeshElement::GeomFilter                                 TFilter;
+    typedef SMDS_SetIterator < TValue, TSetIterator, TAccesor, TFilter > TElemSetIter;
+    SMDS_ElemIteratorPtr elemIt = SMDS_ElemIteratorPtr
+      ( new TElemSetIter( theHexas.begin(),
+                          theHexas.end(),
+                          SMDS_MeshElement::GeomFilter( SMDSGeom_HEXA )));
+
+    SMESH_ElementSearcher* searcher =
+      SMESH_MeshAlgos::GetElementSearcher( *myMesh->GetMeshDS(), elemIt );
+
+    startHex = searcher->FindClosestTo( theFacetNormal.Location(), SMDSAbs_Volume );
+
+    delete searcher;
+
+    if ( !startHex )
+      throw SALOME_Exception( THIS_METHOD "startHex not found");
+  }
+
+  // Select a facet of startHex by theFacetNormal
+
+  SMDS_VolumeTool vTool( startHex );
+  double norm[3], dot, maxDot = 0;
+  int facetID = -1;
+  for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
+    if ( vTool.GetFaceNormal( iF, norm[0], norm[1], norm[2] ))
+    {
+      dot = Abs( theFacetNormal.Direction().Dot( gp_Dir( norm[0], norm[1], norm[2] )));
+      if ( dot > maxDot )
+      {
+        facetID = iF;
+        maxDot = dot;
+      }
+    }
+  if ( facetID < 0 )
+    throw SALOME_Exception( THIS_METHOD "facet of startHex not found");
+
+  // Fill theFacets starting from facetID of startHex
+
+  // facets used for seach of volumes adjacent to already treated ones
+  typedef pair< TFacetOfElem::iterator, int > TElemFacets;
+  typedef map< TVolumeFaceKey, TElemFacets  > TFacetMap;
+  TFacetMap facetsToCheck;
+
+  set<const SMDS_MeshNode*> facetNodes;
+  const SMDS_MeshElement*   curHex;
+
+  const bool allHex = ( theHexas.size() == myMesh->NbHexas() );
+
+  while ( startHex )
+  {
+    // move in two directions from startHex via facetID
+    for ( int is2nd = 0; is2nd < 2; ++is2nd )
+    {
+      curHex       = startHex;
+      int curFacet = facetID;
+      if ( is2nd ) // do not treat startHex twice
+      {
+        vTool.Set( curHex );
+        if ( vTool.IsFreeFace( curFacet, &curHex ))
+        {
+          curHex = 0;
+        }
+        else
+        {
+          vTool.GetFaceNodes( curFacet, facetNodes );
+          vTool.Set( curHex );
+          curFacet = vTool.GetFaceIndex( facetNodes );
+        }
+      }
+      while ( curHex )
+      {
+        // store a facet to split
+        if ( curHex->GetGeomType() != SMDSGeom_HEXA )
+        {
+          theFacets.insert( make_pair( curHex, -1 ));
+          break;
+        }
+        if ( !allHex && !theHexas.count( curHex ))
+          break;
+
+        pair< TFacetOfElem::iterator, bool > facetIt2isNew =
+          theFacets.insert( make_pair( curHex, curFacet ));
+        if ( !facetIt2isNew.second )
+          break;
+
+        // remember not-to-split facets in facetsToCheck
+        int oppFacet = vTool.GetOppFaceIndexOfHex( curFacet );
+        for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
+        {
+          if ( iF == curFacet && iF == oppFacet )
+            continue;
+          TVolumeFaceKey facetKey ( vTool, iF );
+          TElemFacets    elemFacet( facetIt2isNew.first, iF );
+          pair< TFacetMap::iterator, bool > it2isnew =
+            facetsToCheck.insert( make_pair( facetKey, elemFacet ));
+          if ( !it2isnew.second )
+            facetsToCheck.erase( it2isnew.first ); // adjacent hex already checked
+        }
+        // pass to a volume adjacent via oppFacet
+        if ( vTool.IsFreeFace( oppFacet, &curHex ))
+        {
+          curHex = 0;
+        }
+        else
+        {
+          // get a new curFacet
+          vTool.GetFaceNodes( oppFacet, facetNodes );
+          vTool.Set( curHex );
+          curFacet = vTool.GetFaceIndex( facetNodes, /*hint=*/curFacet );
+        }
+      }
+    } // move in two directions from startHex via facetID
+
+    // Find a new startHex by facetsToCheck
+
+    startHex = 0;
+    facetID  = -1;
+    TFacetMap::iterator fIt = facetsToCheck.begin();
+    while ( !startHex && fIt != facetsToCheck.end() )
+    {
+      const TElemFacets&  elemFacets = fIt->second;
+      const SMDS_MeshElement*    hex = elemFacets.first->first;
+      int                 splitFacet = elemFacets.first->second;
+      int               lateralFacet = elemFacets.second;
+      facetsToCheck.erase( fIt );
+      fIt = facetsToCheck.begin();
+
+      vTool.Set( hex );
+      if ( vTool.IsFreeFace( lateralFacet, &curHex ) || 
+           curHex->GetGeomType() != SMDSGeom_HEXA )
+        continue;
+      if ( !allHex && !theHexas.count( curHex ))
+        continue;
+
+      startHex = curHex;
+
+      // find a facet of startHex to split 
+
+      set<const SMDS_MeshNode*> lateralNodes;
+      vTool.GetFaceNodes( lateralFacet, lateralNodes );
+      vTool.GetFaceNodes( splitFacet,   facetNodes );
+      int oppLateralFacet = vTool.GetOppFaceIndexOfHex( lateralFacet );
+      vTool.Set( startHex );
+      lateralFacet = vTool.GetFaceIndex( lateralNodes, oppLateralFacet );
+
+      // look for a facet of startHex having common nodes with facetNodes
+      // but not lateralFacet
+      for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
+      {
+        if ( iF == lateralFacet )
+          continue;
+        int nbCommonNodes = 0;
+        const SMDS_MeshNode** nn = vTool.GetFaceNodes( iF );
+        for ( int iN = 0, nbN = vTool.NbFaceNodes( iF ); iN < nbN; ++iN )
+          nbCommonNodes += facetNodes.count( nn[ iN ]);
+
+        if ( nbCommonNodes >= 2 )
+        {
+          facetID = iF;
+          break;
+        }
+      }
+      if ( facetID < 0 )
+        throw SALOME_Exception( THIS_METHOD "facet of a new startHex not found");
+    }
+  } //   while ( startHex )
+}
+
 //=======================================================================
 //function : AddToSameGroups
 //purpose  : add elemToAdd to the groups the elemInGroups belongs to
@@ -3200,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 );
@@ -3275,8 +3768,11 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
   // smooth elements on each TopoDS_Face separately
   // ===============================================
 
-  set< int >::reverse_iterator fId = faceIdSet.rbegin(); // treate 0 fId at the end
-  for ( ; fId != faceIdSet.rend(); ++fId ) {
+  SMESH_MesherHelper helper( *GetMesh() );
+
+  set< int >::reverse_iterator fId = faceIdSet.rbegin(); // treat 0 fId at the end
+  for ( ; fId != faceIdSet.rend(); ++fId )
+  {
     // get face surface and submesh
     Handle(Geom_Surface) surface;
     SMESHDS_SubMesh* faceSubMesh = 0;
@@ -3284,7 +3780,8 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
     double fToler2 = 0, f,l;
     double u1 = 0, u2 = 0, v1 = 0, v2 = 0;
     bool isUPeriodic = false, isVPeriodic = false;
-    if ( *fId ) {
+    if ( *fId )
+    {
       face = TopoDS::Face( aMesh->IndexToShape( *fId ));
       surface = BRep_Tool::Surface( face );
       faceSubMesh = aMesh->MeshElements( *fId );
@@ -3297,6 +3794,7 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
       if ( isVPeriodic )
         surface->VPeriod();
       surface->Bounds( u1, u2, v1, v2 );
+      helper.SetSubShape( face );
     }
     // ---------------------------------------------------------
     // for elements on a face, find movable and fixed nodes and
@@ -3318,7 +3816,8 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
     int nbElemOnFace = 0;
     itElem = theElems.begin();
     // loop on not yet smoothed elements: look for elems on a face
-    while ( itElem != theElems.end() ) {
+    while ( itElem != theElems.end() )
+    {
       if ( faceSubMesh && nbElemOnFace == faceSubMesh->NbElements() )
         break; // all elements found
 
@@ -3374,12 +3873,15 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
 
       // get nodes to check UV
       list< const SMDS_MeshNode* > uvCheckNodes;
+      const SMDS_MeshNode* nodeInFace = 0;
       itN = elem->nodesIterator();
       nn = 0; nbn =  elem->NbNodes();
       if(elem->IsQuadratic())
         nbn = nbn/2;
       while ( nn++ < nbn ) {
         node = static_cast<const SMDS_MeshNode*>( itN->next() );
+        if ( node->GetPosition()->GetDim() == 2 )
+          nodeInFace = node;
         if ( uvMap.find( node ) == uvMap.end() )
           uvCheckNodes.push_back( node );
         // add nodes of elems sharing node
@@ -3405,41 +3907,21 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
         const SMDS_PositionPtr& pos = node->GetPosition();
         posType = pos ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
         // get existing UV
-        switch ( posType ) {
-        case SMDS_TOP_FACE: {
-          SMDS_FacePosition* fPos = ( SMDS_FacePosition* ) pos;
-          uv.SetCoord( fPos->GetUParameter(), fPos->GetVParameter() );
-          break;
-        }
-        case SMDS_TOP_EDGE: {
-          TopoDS_Shape S = aMesh->IndexToShape( node->getshapeId() );
-          Handle(Geom2d_Curve) pcurve;
-          if ( !S.IsNull() && S.ShapeType() == TopAbs_EDGE )
-            pcurve = BRep_Tool::CurveOnSurface( TopoDS::Edge( S ), face, f,l );
-          if ( !pcurve.IsNull() ) {
-            double u = (( SMDS_EdgePosition* ) pos )->GetUParameter();
-            uv = pcurve->Value( u ).XY();
-          }
-          break;
-        }
-        case SMDS_TOP_VERTEX: {
-          TopoDS_Shape S = aMesh->IndexToShape( node->getshapeId() );
-          if ( !S.IsNull() && S.ShapeType() == TopAbs_VERTEX )
-            uv = BRep_Tool::Parameters( TopoDS::Vertex( S ), face ).XY();
-          break;
-        }
-        default:;
-        }
-        // check existing UV
-        bool project = true;
-        gp_Pnt pNode ( node->X(), node->Y(), node->Z() );
-        double dist1 = DBL_MAX, dist2 = 0;
-        if ( posType != SMDS_TOP_3DSPACE ) {
-          dist1 = pNode.SquareDistance( surface->Value( uv.X(), uv.Y() ));
-          project = dist1 > fToler2;
-        }
+        if ( pos )
+        {
+          bool toCheck = true;
+          uv = helper.GetNodeUV( face, node, nodeInFace, &toCheck );
+        }
+        // compute not existing UV
+        bool project = ( posType == SMDS_TOP_3DSPACE );
+        // double dist1 = DBL_MAX, dist2 = 0;
+        // if ( posType != SMDS_TOP_3DSPACE ) {
+        //   dist1 = pNode.SquareDistance( surface->Value( uv.X(), uv.Y() ));
+        //   project = dist1 > fToler2;
+        // }
         if ( project ) { // compute new UV
           gp_XY newUV;
+          gp_Pnt pNode = SMESH_TNodeXYZ( node );
           if ( !getClosestUV( projector, pNode, newUV )) {
             MESSAGE("Node Projection Failed " << node);
           }
@@ -3449,9 +3931,9 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
             if ( isVPeriodic )
               newUV.SetY( ElCLib::InPeriod( newUV.Y(), v1, v2 ));
             // check new UV
-            if ( posType != SMDS_TOP_3DSPACE )
-              dist2 = pNode.SquareDistance( surface->Value( newUV.X(), newUV.Y() ));
-            if ( dist2 < dist1 )
+            // if ( posType != SMDS_TOP_3DSPACE )
+            //   dist2 = pNode.SquareDistance( surface->Value( newUV.X(), newUV.Y() ));
+            // if ( dist2 < dist1 )
               uv = newUV;
           }
         }
@@ -3519,9 +4001,8 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
         uv2 = pcurve->Value( f );
         int iPar = Abs( uv1.X() - uv2.X() ) > Abs( uv1.Y() - uv2.Y() ) ? 1 : 2;
         // assure uv1 < uv2
-        if ( uv1.Coord( iPar ) > uv2.Coord( iPar )) {
-          gp_Pnt2d tmp = uv1; uv1 = uv2; uv2 = tmp;
-        }
+        if ( uv1.Coord( iPar ) > uv2.Coord( iPar ))
+          std::swap( uv1, uv2 );
         // get nodes on seam and its vertices
         list< const SMDS_MeshNode* > seamNodes;
         SMDS_NodeIteratorPtr nSeamIt = sm->GetNodes();
@@ -3571,12 +4052,14 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
                   setMovableNodes.find( n ) == setMovableNodes.end() )
                 continue;
               // add only nodes being closer to uv2 than to uv1
-              gp_Pnt pMid (0.5 * ( n->X() + nSeam->X() ),
-                           0.5 * ( n->Y() + nSeam->Y() ),
-                           0.5 * ( n->Z() + nSeam->Z() ));
-              gp_XY uv;
-              getClosestUV( projector, pMid, uv );
-              if ( uv.Coord( iPar ) > uvMap[ n ]->Coord( iPar ) ) {
+              // gp_Pnt pMid (0.5 * ( n->X() + nSeam->X() ),
+              //              0.5 * ( n->Y() + nSeam->Y() ),
+              //              0.5 * ( n->Z() + nSeam->Z() ));
+              // gp_XY uv;
+              // getClosestUV( projector, pMid, uv );
+              double x = uvMap[ n ]->Coord( iPar );
+              if ( Abs( uv1.Coord( iPar ) - x ) >
+                   Abs( uv2.Coord( iPar ) - x )) {
                 nodesNearSeam.insert( n );
                 nbUseMap2++;
               }
@@ -3685,8 +4168,6 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
     // move medium nodes of quadratic elements
     if ( isQuadratic )
     {
-      SMESH_MesherHelper helper( *GetMesh() );
-      helper.SetSubShape( face );
       vector<const SMDS_MeshNode*> nodes;
       bool checkUV;
       list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
@@ -3723,26 +4204,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;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Assure that theElemSets[0] holds elements, not nodes
+   */
+  //================================================================================
 
-  return faceNorm * extrDir < 0.0;
+  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] );
+    }
+  }
 }
 
 //=======================================================================
@@ -3842,6 +4346,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 ) {
@@ -3947,11 +4475,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 )
@@ -4117,8 +4645,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
     {
@@ -4172,7 +4700,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 );
@@ -4184,7 +4713,7 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
     for ( iNode = 0; iNode < nbNodes; iNode++ )
       prevNod[ iNode ] = nextNod[ iNode ];
 
-  } // for steps
+  } // loop on steps
 }
 
 //=======================================================================
@@ -4594,177 +5123,460 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
           srcElements.Append( elem );
       }
     }
-  } // loop on swept elements
+  } // loop on swept elements
+}
+
+//=======================================================================
+//function : RotationSweep
+//purpose  :
+//=======================================================================
+
+SMESH_MeshEditor::PGroupIDs
+SMESH_MeshEditor::RotationSweep(TIDSortedElemSet   theElemSets[2],
+                                const gp_Ax1&      theAxis,
+                                const double       theAngle,
+                                const int          theNbSteps,
+                                const double       theTol,
+                                const bool         theMakeGroups,
+                                const bool         theMakeWalls)
+{
+  myLastCreatedElems.Clear();
+  myLastCreatedNodes.Clear();
+
+  // source elements for each generated one
+  SMESH_SequenceOfElemPtr srcElems, srcNodes;
+
+  MESSAGE( "RotationSweep()");
+  gp_Trsf aTrsf;
+  aTrsf.SetRotation( theAxis, theAngle );
+  gp_Trsf aTrsf2;
+  aTrsf2.SetRotation( theAxis, theAngle/2. );
+
+  gp_Lin aLine( theAxis );
+  double aSqTol = theTol * theTol;
+
+  SMESHDS_Mesh* aMesh = GetMeshDS();
+
+  TNodeOfNodeListMap mapNewNodes;
+  TElemOfVecOfNnlmiMap mapElemNewNodes;
+  TTElemOfElemListMap newElemsMap;
+
+  const bool isQuadraticMesh = bool( myMesh->NbEdges(ORDER_QUADRATIC) +
+                                     myMesh->NbFaces(ORDER_QUADRATIC) +
+                                     myMesh->NbVolumes(ORDER_QUADRATIC) );
+  // loop on theElemSets
+  setElemsFirst( theElemSets );
+  TIDSortedElemSet::iterator itElem;
+  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() );
+
+        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 )
+          {
+            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] );
+                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 );
+            }
+            else {
+              listNewNodes.push_back( newNode );
+              // if ( needMediumNodes )
+              //   listNewNodes.push_back( newNode );
+            }
+          }
+        }
+        newNodesItVec.push_back( nIt );
+      }
+      // make new elements
+      sweepElement( elem, newNodesItVec, newElemsMap[elem], theNbSteps, srcElems );
+    }
+  }
+
+  if ( theMakeWalls )
+    makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElemSets[0], theNbSteps, srcElems );
+
+  PGroupIDs newGroupIDs;
+  if ( theMakeGroups )
+    newGroupIDs = generateGroups( srcNodes, srcElems, "rotated");
+
+  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 : ExtrusParam
+//purpose  : steps are given explicitly
+//=======================================================================
+
+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 )
+{
+  if (( theFlags & EXTRUSION_FLAG_SEW ) &&
+      ( theTolerance > 0 ))
+  {
+    myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByDirAndSew;
+  }
+  else
+  {
+    myMakeNodesFun = & SMESH_MeshEditor::ExtrusParam::makeNodesByDir;
+  }
+}
+
+//=======================================================================
+//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;
+  }
+}
+
+//=======================================================================
+//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();
+  }
+  return res;
 }
 
 //=======================================================================
-//function : RotationSweep
-//purpose  :
+//function : ExtrusParam::makeNodesByDir
+//purpose  : create nodes for standard extrusion
 //=======================================================================
 
-SMESH_MeshEditor::PGroupIDs
-SMESH_MeshEditor::RotationSweep(TIDSortedElemSet & theElems,
-                                const gp_Ax1&      theAxis,
-                                const double       theAngle,
-                                const int          theNbSteps,
-                                const double       theTol,
-                                const bool         theMakeGroups,
-                                const bool         theMakeWalls)
+int SMESH_MeshEditor::ExtrusParam::
+makeNodesByDir( SMESHDS_Mesh*                     mesh,
+                const SMDS_MeshNode*              srcNode,
+                std::list<const SMDS_MeshNode*> & newNodes,
+                const bool                        makeMediumNodes)
 {
-  myLastCreatedElems.Clear();
-  myLastCreatedNodes.Clear();
-
-  // source elements for each generated one
-  SMESH_SequenceOfElemPtr srcElems, srcNodes;
-
-  MESSAGE( "RotationSweep()");
-  gp_Trsf aTrsf;
-  aTrsf.SetRotation( theAxis, theAngle );
-  gp_Trsf aTrsf2;
-  aTrsf2.SetRotation( theAxis, theAngle/2. );
-
-  gp_Lin aLine( theAxis );
-  double aSqTol = theTol * theTol;
-
-  SMESHDS_Mesh* aMesh = GetMeshDS();
+  gp_XYZ p = SMESH_TNodeXYZ( srcNode );
 
-  TNodeOfNodeListMap mapNewNodes;
-  TElemOfVecOfNnlmiMap mapElemNewNodes;
-  TTElemOfElemListMap newElemsMap;
-
-  const bool isQuadraticMesh = bool( myMesh->NbEdges(ORDER_QUADRATIC) +
-                                     myMesh->NbFaces(ORDER_QUADRATIC) +
-                                     myMesh->NbVolumes(ORDER_QUADRATIC) );
-  // loop on theElems
-  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() );
+  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;
+}
 
-    // 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() );
+//=======================================================================
+//function : ExtrusParam::makeNodesByDirAndSew
+//purpose  : create nodes for standard extrusion with sewing
+//=======================================================================
 
-      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 );
+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 );
 
-      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 )
+  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 )
         {
-          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;
-          }
+          node = myNodes.Value(i);
+          break;
         }
-
-        // 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 );
-          }
-          else {
-            listNewNodes.push_back( newNode );
-            // if ( needMediumNodes )
-            //   listNewNodes.push_back( newNode );
-          }
+      }
+    }
+    else {
+      SMDS_NodeIteratorPtr itn = mesh->nodesIterator();
+      while(itn->more()) {
+        SMESH_TNodeXYZ P2( itn->next() );
+        if (( P1 - P2 ).SquareModulus() < myTolerance * myTolerance )
+        {
+          node = P2._node;
+          break;
         }
       }
-      newNodesItVec.push_back( nIt );
     }
-    // make new elements
-    sweepElement( elem, newNodesItVec, newElemsMap[elem], theNbSteps, srcElems );
-  }
 
-  if ( theMakeWalls )
-    makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElems, theNbSteps, srcElems );
+    if ( !node )
+      node = mesh->AddNode( P1.X(), P1.Y(), P1.Z() );
 
-  PGroupIDs newGroupIDs;
-  if ( theMakeGroups )
-    newGroupIDs = generateGroups( srcNodes, srcElems, "rotated");
+    newNodes.push_back( node );
 
-  return newGroupIDs;
-}
+  } // loop on steps
 
+  return nbNodes;
+}
 
 //=======================================================================
-//function : CreateNode
-//purpose  :
+//function : ExtrusParam::makeNodesByNormal2D
+//purpose  : create nodes for extrusion using normals of faces
 //=======================================================================
-const SMDS_MeshNode* SMESH_MeshEditor::CreateNode(const double x,
-                                                  const double y,
-                                                  const double z,
-                                                  const double tolnode,
-                                                  SMESH_SequenceOfNode& aNodes)
+
+int SMESH_MeshEditor::ExtrusParam::
+makeNodesByNormal2D( SMESHDS_Mesh*                     mesh,
+                     const SMDS_MeshNode*              srcNode,
+                     std::list<const SMDS_MeshNode*> & newNodes,
+                     const bool                        makeMediumNodes)
 {
-  // myLastCreatedElems.Clear();
-  // myLastCreatedNodes.Clear();
+  const bool alongAvgNorm = ( myFlags & EXTRUSION_FLAG_BY_AVG_NORMAL );
 
-  gp_Pnt P1(x,y,z);
-  SMESHDS_Mesh * aMesh = myMesh->GetMeshDS();
+  gp_XYZ p = SMESH_TNodeXYZ( srcNode );
 
-  // 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);
+  // 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 );
+      }
     }
   }
-  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;
-    }
+
+  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 );
   }
 
-  // create new node and return it
-  const SMDS_MeshNode* NewNode = aMesh->AddNode(x,y,z);
-  //myLastCreatedNodes.Append(NewNode);
-  return NewNode;
+  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
@@ -4772,24 +5584,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 );
 }
 
 
@@ -4799,12 +5602,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();
@@ -4814,7 +5614,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;
@@ -4826,91 +5628,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;
@@ -4921,7 +5707,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,
@@ -4950,7 +5736,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;
@@ -5093,7 +5879,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,
@@ -5121,7 +5907,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;
@@ -5162,7 +5948,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
     }
 
     conn = nbEdgeConnectivity(theN1);
-    if(conn > 2)
+    if( conn != 1 )
       return EXTR_PATH_NOT_EDGE;
 
     aItE = theN1->GetInverseElementIterator();
@@ -5442,7 +6228,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,
@@ -5480,15 +6266,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();
@@ -5498,144 +6288,151 @@ 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;
+      if ( !elem )
+        continue;
+      // SMDSAbs_ElementType aTypeE = elem->GetType();
+      // if ( 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");
@@ -6011,7 +6808,7 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
 
   if ( ( theMakeGroups && theCopy ) ||
        ( theMakeGroups && theTargetMesh ) )
-    newGroupIDs = generateGroups( srcNodes, srcElems, groupPostfix, theTargetMesh );
+    newGroupIDs = generateGroups( srcNodes, srcElems, groupPostfix, theTargetMesh, false );
 
   return newGroupIDs;
 }
@@ -6019,9 +6816,11 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
 //=======================================================================
 /*!
  * \brief Create groups of elements made during transformation
- * \param nodeGens - nodes making corresponding myLastCreatedNodes
- * \param elemGens - elements making corresponding myLastCreatedElems
- * \param postfix - to append to names of new groups
+ *  \param nodeGens - nodes making corresponding myLastCreatedNodes
+ *  \param elemGens - elements making corresponding myLastCreatedElems
+ *  \param postfix - to append to names of new groups
+ *  \param targetMesh - mesh to create groups in
+ *  \param topPresent - is there "top" elements that are created by sweeping
  */
 //=======================================================================
 
@@ -6029,14 +6828,17 @@ SMESH_MeshEditor::PGroupIDs
 SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
                                  const SMESH_SequenceOfElemPtr& elemGens,
                                  const std::string&             postfix,
-                                 SMESH_Mesh*                    targetMesh)
+                                 SMESH_Mesh*                    targetMesh,
+                                 const bool                     topPresent)
 {
   PGroupIDs newGroupIDs( new list<int> );
   SMESH_Mesh* mesh = targetMesh ? targetMesh : GetMesh();
 
   // Sort existing groups by types and collect their names
 
-  // to store an old group and a generated new ones
+  // containers to store an old group and generated new ones;
+  // 1st new group is for result elems of different type than a source one;
+  // 2nd new group is for same type result elems ("top" group at extrusion)
   using boost::tuple;
   using boost::make_tuple;
   typedef tuple< SMESHDS_GroupBase*, SMESHDS_Group*, SMESHDS_Group* > TOldNewGroup;
@@ -6066,6 +6868,7 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
 
   // Loop on nodes and elements to add them in new groups
 
+  vector< const SMDS_MeshElement* > resultElems;
   for ( int isNodes = 0; isNodes < 2; ++isNodes )
   {
     const SMESH_SequenceOfElemPtr& gens  = isNodes ? nodeGens : elemGens;
@@ -6088,7 +6891,7 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
         continue;
       }
       // collect all elements made by the iElem-th sourceElem
-      list< const SMDS_MeshElement* > resultElems;
+      resultElems.clear();
       if ( const SMDS_MeshElement* resElem = elems( iElem ))
         if ( resElem != sourceElem )
           resultElems.push_back( resElem );
@@ -6097,25 +6900,23 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
           if ( resElem != sourceElem )
             resultElems.push_back( resElem );
 
-      // there must be a top element
       const SMDS_MeshElement* topElem = 0;
-      if ( isNodes )
+      if ( isNodes ) // there must be a top element
       {
         topElem = resultElems.back();
         resultElems.pop_back();
       }
       else
       {
-        list< const SMDS_MeshElement* >::reverse_iterator resElemIt = resultElems.rbegin();
+        vector< const SMDS_MeshElement* >::reverse_iterator resElemIt = resultElems.rbegin();
         for ( ; resElemIt != resultElems.rend() ; ++resElemIt )
           if ( (*resElemIt)->GetType() == sourceElem->GetType() )
           {
             topElem = *resElemIt;
-            resultElems.erase( --(resElemIt.base()) ); // erase *resElemIt
+            *resElemIt = 0; // erase *resElemIt
             break;
           }
       }
-
       // add resultElems to groups originted from ones the sourceElem belongs to
       list< TOldNewGroup >::iterator gOldNew, gLast = groupsOldNew.end();
       for ( gOldNew = groupsOldNew.begin(); gOldNew != gLast; ++gOldNew )
@@ -6125,16 +6926,17 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
         {
           // fill in a new group
           SMDS_MeshGroup & newGroup = gOldNew->get<1>()->SMDSGroup();
-          list< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt;
+          vector< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt;
           for ( resElemIt = resultElems.begin(); resElemIt != resLast; ++resElemIt )
-            newGroup.Add( *resElemIt );
+            if ( *resElemIt )
+              newGroup.Add( *resElemIt );
 
           // fill a "top" group
           if ( topElem )
           {
             SMDS_MeshGroup & newTopGroup = gOldNew->get<2>()->SMDSGroup();
             newTopGroup.Add( topElem );
-          }
+         }
         }
       }
     } // loop on created elements
@@ -6148,7 +6950,6 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
     SMESHDS_GroupBase* oldGroupDS =   orderedOldNewGroups[i]->get<0>();
     SMESHDS_Group*   newGroups[2] = { orderedOldNewGroups[i]->get<1>(),
                                       orderedOldNewGroups[i]->get<2>() };
-    const int nbNewGroups = !newGroups[0]->IsEmpty() + !newGroups[1]->IsEmpty();
     for ( int is2nd = 0; is2nd < 2; ++is2nd )
     {
       SMESHDS_Group* newGroupDS = newGroups[ is2nd ];
@@ -6162,11 +6963,21 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
         newGroupDS->SetType( newGroupDS->GetElements()->next()->GetType() );
 
         // make a name
-        const bool isTop = ( nbNewGroups == 2 &&
+        const bool isTop = ( topPresent &&
                              newGroupDS->GetType() == oldGroupDS->GetType() &&
                              is2nd );
 
         string name = oldGroupDS->GetStoreName();
+        { // remove trailing whitespaces (issue 22599)
+          size_t size = name.size();
+          while ( size > 1 && isspace( name[ size-1 ]))
+            --size;
+          if ( size != name.size() )
+          {
+            name.resize( size );
+            oldGroupDS->SetStoreName( name.c_str() );
+          }
+        }
         if ( !targetMesh ) {
           string suffix = ( isTop ? "top": postfix.c_str() );
           name += "_";
@@ -8334,6 +9145,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT
   aHelper.SetIsQuadratic( true );
   aHelper.SetIsBiQuadratic( theToBiQuad );
   aHelper.SetElementsOnShape(true);
+  aHelper.ToFixNodeParameters( true );
 
   // convert elements assigned to sub-meshes
   int nbCheckedElems = 0;
@@ -9558,7 +10370,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
@@ -9717,6 +10529,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 );
       }
@@ -9729,10 +10542,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;
@@ -9743,8 +10554,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
 */
 //================================================================================
@@ -9780,6 +10591,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 );
     }
@@ -9881,15 +10693,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
     {
@@ -10146,7 +10955,12 @@ double SMESH_MeshEditor::OrientedAngle(const gp_Pnt& p0, const gp_Pnt& p1, const
   gp_Vec v2(p0, g2);
   gp_Vec n1 = vref.Crossed(v1);
   gp_Vec n2 = vref.Crossed(v2);
-  return n2.AngleWithRef(n1, vref);
+  try {
+    return n2.AngleWithRef(n1, vref);
+  }
+  catch ( Standard_Failure ) {
+  }
+  return Max( v1.Magnitude(), v2.Magnitude() );
 }
 
 /*!
@@ -10159,13 +10973,16 @@ double SMESH_MeshEditor::OrientedAngle(const gp_Pnt& p0, const gp_Pnt& p1, const
  * If there is no shared faces between the group #n and the group #p in the list, the group j_n_p is not created.
  * All the flat elements are gathered into the group named "joints3D" (or "joints2D" in 2D situation).
  * The flat element of the multiple junctions between the simple junction are stored in a group named "jointsMultiples".
- * @param theElems - list of groups of volumes, where a group of volume is a set of
- * SMDS_MeshElements sorted by Id.
- * @param createJointElems - if TRUE, create the elements
- * @return TRUE if operation has been completed successfully, FALSE otherwise
+ * \param theElems - list of groups of volumes, where a group of volume is a set of
+ *        SMDS_MeshElements sorted by Id.
+ * \param createJointElems - if TRUE, create the elements
+ * \param onAllBoundaries - if TRUE, the nodes and elements are also created on
+ *        the boundary between \a theDomains and the rest mesh
+ * \return TRUE if operation has been completed successfully, FALSE otherwise
  */
 bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSortedElemSet>& theElems,
-                                                     bool createJointElems)
+                                                     bool                                 createJointElems,
+                                                     bool                                 onAllBoundaries)
 {
   MESSAGE("----------------------------------------------");
   MESSAGE("SMESH_MeshEditor::doubleNodesOnGroupBoundaries");
@@ -10194,15 +11011,20 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
 
   MESSAGE(".. Number of domains :"<<theElems.size());
 
+  TIDSortedElemSet theRestDomElems;
+  const int iRestDom  = -1;
+  const int idom0     = onAllBoundaries ? iRestDom : 0;
+  const int nbDomains = theElems.size();
+
   // Check if the domains do not share an element
-  for (int idom = 0; idom < theElems.size()-1; idom++)
+  for (int idom = 0; idom < nbDomains-1; idom++)
     {
 //       MESSAGE("... Check of domain #" << idom);
       const TIDSortedElemSet& domain = theElems[idom];
       TIDSortedElemSet::const_iterator elemItr = domain.begin();
       for (; elemItr != domain.end(); ++elemItr)
         {
-          SMDS_MeshElement* anElem = (SMDS_MeshElement*) *elemItr;
+          const SMDS_MeshElement* anElem = *elemItr;
           int idombisdeb = idom + 1 ;
           for (int idombis = idombisdeb; idombis < theElems.size(); idombis++) // check if the element belongs to a domain further in the list
           {
@@ -10218,7 +11040,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
         }
     }
 
-  for (int idom = 0; idom < theElems.size(); idom++)
+  for (int idom = 0; idom < nbDomains; idom++)
     {
 
       // --- build a map (face to duplicate --> volume to modify)
@@ -10231,7 +11053,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       TIDSortedElemSet::const_iterator elemItr = domain.begin();
       for (; elemItr != domain.end(); ++elemItr)
         {
-          SMDS_MeshElement* anElem = (SMDS_MeshElement*) *elemItr;
+          const SMDS_MeshElement* anElem = *elemItr;
           if (!anElem)
             continue;
           int vtkId = anElem->getVtkId();
@@ -10244,26 +11066,30 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
             {
               int smdsId = meshDS->fromVtkToSmds(neighborsVtkIds[n]);
               const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
-              if (! domain.count(elem)) // neighbor is in another domain : face is shared
+              if (elem && ! domain.count(elem)) // neighbor is in another domain : face is shared
                 {
                   bool ok = false ;
-                  for (int idombis = 0; idombis < theElems.size(); idombis++) // check if the neighbor belongs to another domain of the list
+                  for (int idombis = 0; idombis < theElems.size() && !ok; idombis++) // check if the neighbor belongs to another domain of the list
                   {
                     // MESSAGE("Domain " << idombis);
                     const TIDSortedElemSet& domainbis = theElems[idombis];
                     if ( domainbis.count(elem)) ok = true ; // neighbor is in a correct domain : face is kept
                   }
-                  if ( ok ) // the characteristics of the face is stored
+                  if ( ok || onAllBoundaries ) // the characteristics of the face is stored
                   {
                     DownIdType face(downIds[n], downTypes[n]);
-                    if (!faceDomains.count(face))
-                      faceDomains[face] = emptyMap; // create an empty entry for face
                     if (!faceDomains[face].count(idom))
                       {
                         faceDomains[face][idom] = vtkId; // volume associated to face in this domain
                         celldom[vtkId] = idom;
                         //MESSAGE("       cell with a border " << vtkId << " domain " << idom);
                       }
+                    if ( !ok )
+                    {
+                      theRestDomElems.insert( elem );
+                      faceDomains[face][iRestDom] = neighborsVtkIds[n];
+                      celldom[neighborsVtkIds[n]] = iRestDom;
+                    }
                   }
                 }
             }
@@ -10277,14 +11103,14 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   //     explore the nodes of the face and see if they belong to a cell in the domain,
   //     which has only a node or an edge on the border (not a shared face)
 
-  for (int idomain = 0; idomain < theElems.size(); idomain++)
+  for (int idomain = idom0; idomain < nbDomains; idomain++)
     {
       //MESSAGE("Domain " << idomain);
-      const TIDSortedElemSet& domain = theElems[idomain];
+      const TIDSortedElemSet& domain = (idomain == iRestDom) ? theRestDomElems : theElems[idomain];
       itface = faceDomains.begin();
       for (; itface != faceDomains.end(); ++itface)
         {
-          std::map<int, int> domvol = itface->second;
+          const std::map<int, int>& domvol = itface->second;
           if (!domvol.count(idomain))
             continue;
           DownIdType face = itface->first;
@@ -10313,8 +11139,6 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
                                 //no cells created after BuildDownWardConnectivity
                     }
                   DownIdType aCell(downId, vtkType);
-                  if (!cellDomains.count(aCell))
-                    cellDomains[aCell] = emptyMap; // create an empty entry for cell
                   cellDomains[aCell][idomain] = vtkId;
                   celldom[vtkId] = idomain;
                   //MESSAGE("       cell " << vtkId << " domain " << idomain);
@@ -10336,12 +11160,12 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   std::map<int, std::vector<int> > mutipleNodesToFace; // nodes multi domains with domain order to transform in Face (junction between 3 or more 2D domains)
 
   MESSAGE(".. Duplication of the nodes");
-  for (int idomain = 0; idomain < theElems.size(); idomain++)
+  for (int idomain = idom0; idomain < nbDomains; idomain++)
     {
       itface = faceDomains.begin();
       for (; itface != faceDomains.end(); ++itface)
         {
-          std::map<int, int> domvol = itface->second;
+          const std::map<int, int>& domvol = itface->second;
           if (!domvol.count(idomain))
             continue;
           DownIdType face = itface->first;
@@ -10353,15 +11177,12 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
           for (; itn != oldNodes.end(); ++itn)
             {
               int oldId = *itn;
-              //MESSAGE("-+-+-a node " << oldId);
-              if (!nodeDomains.count(oldId))
-                nodeDomains[oldId] = emptyMap; // create an empty entry for node
               if (nodeDomains[oldId].empty())
                 {
                   nodeDomains[oldId][idomain] = oldId; // keep the old node in the first domain
                   //MESSAGE("-+-+-b     oldNode " << oldId << " domain " << idomain);
                 }
-              std::map<int, int>::iterator itdom = domvol.begin();
+              std::map<int, int>::const_iterator itdom = domvol.begin();
               for (; itdom != domvol.end(); ++itdom)
                 {
                   int idom = itdom->first;
@@ -10389,6 +11210,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());
@@ -10399,7 +11221,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
     }
 
   MESSAGE(".. Creation of elements");
-  for (int idomain = 0; idomain < theElems.size(); idomain++)
+  for (int idomain = idom0; idomain < nbDomains; idomain++)
     {
       itface = faceDomains.begin();
       for (; itface != faceDomains.end(); ++itface)
@@ -10468,11 +11290,12 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
                               for (int id=0; id < doms.size(); id++)
                                 {
                                   int idom = doms[id];
+                                  const TIDSortedElemSet& domain = (idom == iRestDom) ? theRestDomElems : theElems[idom];
                                   for (int ivol=0; ivol<nbvol; ivol++)
                                     {
                                       int smdsId = meshDS->fromVtkToSmds(vtkVolIds[ivol]);
                                       SMDS_MeshElement* elem = (SMDS_MeshElement*)meshDS->FindElement(smdsId);
-                                      if (theElems[idom].count(elem))
+                                      if (domain.count(elem))
                                         {
                                           SMDS_VtkVolume* svol = dynamic_cast<SMDS_VtkVolume*>(elem);
                                           domvol[idom] = svol;
@@ -10653,7 +11476,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   feDom.clear();
 
   MESSAGE(".. Modification of elements");
-  for (int idomain = 0; idomain < theElems.size(); idomain++)
+  for (int idomain = idom0; idomain < nbDomains; idomain++)
     {
       std::map<int, std::map<int, int> >::const_iterator itnod = nodeDomains.begin();
       for (; itnod != nodeDomains.end(); ++itnod)
@@ -10726,6 +11549,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();
 
@@ -10793,6 +11624,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
@@ -10809,6 +11641,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
@@ -11753,3 +12586,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:;
+  }
+}