Salome HOME
Merge from V5_1_main 10/12/2010
[modules/smesh.git] / src / SMESH / SMESH_MeshEditor.cxx
index 62de976cfe04d6cc164f27cec0685d2ae1fc5c42..d21d4784c712103ba63f4ad84efb7e5d0b4af120 100644 (file)
@@ -1,4 +1,4 @@
-//  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
+//  Copyright (C) 2007-2010  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
 //
 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
+
 //  SMESH SMESH : idl implementation based on 'SMESH' unit's classes
 // File      : SMESH_MeshEditor.cxx
 // Created   : Mon Apr 12 16:10:22 2004
 // Author    : Edward AGAPOV (eap)
 //
+#define CHRONODEF
 #include "SMESH_MeshEditor.hxx"
 
 #include "SMDS_FaceOfNodes.hxx"
 #include "SMDS_PolyhedralVolumeOfNodes.hxx"
 #include "SMDS_FacePosition.hxx"
 #include "SMDS_SpacePosition.hxx"
-#include "SMDS_QuadraticFaceOfNodes.hxx"
+//#include "SMDS_QuadraticFaceOfNodes.hxx"
 #include "SMDS_MeshGroup.hxx"
+#include "SMDS_LinearEdge.hxx"
+#include "SMDS_Downward.hxx"
+#include "SMDS_SetIterator.hxx"
 
 #include "SMESHDS_Group.hxx"
 #include "SMESHDS_Mesh.hxx"
 
-#include "SMESH_subMesh.hxx"
+#include "SMESH_Algo.hxx"
 #include "SMESH_ControlsDef.hxx"
+#include "SMESH_Group.hxx"
 #include "SMESH_MesherHelper.hxx"
 #include "SMESH_OctreeNode.hxx"
-#include "SMESH_Group.hxx"
+#include "SMESH_subMesh.hxx"
 
 #include "utilities.h"
 
-#include <BRep_Tool.hxx>
+#include <BRepAdaptor_Surface.hxx>
 #include <BRepClass3d_SolidClassifier.hxx>
+#include <BRep_Tool.hxx>
 #include <ElCLib.hxx>
 #include <Extrema_GenExtPS.hxx>
+#include <Extrema_POnCurv.hxx>
 #include <Extrema_POnSurf.hxx>
+#include <GC_MakeSegment.hxx>
 #include <Geom2d_Curve.hxx>
+#include <GeomAPI_ExtremaCurveCurve.hxx>
 #include <GeomAdaptor_Surface.hxx>
 #include <Geom_Curve.hxx>
+#include <Geom_Line.hxx>
 #include <Geom_Surface.hxx>
+#include <IntAna_IntConicQuad.hxx>
+#include <IntAna_Quadric.hxx>
 #include <Precision.hxx>
 #include <TColStd_ListOfInteger.hxx>
 #include <TopAbs_State.hxx>
@@ -74,6 +87,7 @@
 #include <gp_Vec.hxx>
 #include <gp_XY.hxx>
 #include <gp_XYZ.hxx>
+
 #include <math.h>
 
 #include <map>
@@ -88,28 +102,8 @@ using namespace SMESH::Controls;
 
 typedef map<const SMDS_MeshElement*, list<const SMDS_MeshNode*> >    TElemOfNodeListMap;
 typedef map<const SMDS_MeshElement*, list<const SMDS_MeshElement*> > TElemOfElemListMap;
-//typedef map<const SMDS_MeshNode*, vector<const SMDS_MeshNode*> >     TNodeOfNodeVecMap;
-//typedef TNodeOfNodeVecMap::iterator                                  TNodeOfNodeVecMapItr;
-//typedef map<const SMDS_MeshElement*, vector<TNodeOfNodeVecMapItr> >  TElemOfVecOfMapNodesMap;
 
-//=======================================================================
-/*!
- * \brief SMDS_MeshNode -> gp_XYZ convertor
- */
-//=======================================================================
-
-struct TNodeXYZ : public gp_XYZ
-{
-  TNodeXYZ( const SMDS_MeshNode* n ):gp_XYZ( n->X(), n->Y(), n->Z() ) {}
-  double Distance( const SMDS_MeshNode* n )
-  {
-    return gp_Vec( *this, TNodeXYZ( n )).Magnitude();
-  }
-  double SquareDistance( const SMDS_MeshNode* n )
-  {
-    return gp_Vec( *this, TNodeXYZ( n )).SquareMagnitude();
-  }
-};
+typedef SMDS_SetIterator< SMDS_pElement, TIDSortedElemSet::const_iterator> TSetIterator;
 
 //=======================================================================
 //function : SMESH_MeshEditor
@@ -133,99 +127,132 @@ SMESH_MeshEditor::AddElement(const vector<const SMDS_MeshNode*> & node,
                              const bool                           isPoly,
                              const int                            ID)
 {
+  //MESSAGE("AddElement " <<node.size() << " " << type << " " << isPoly << " " << ID);
   SMDS_MeshElement* e = 0;
   int nbnode = node.size();
   SMESHDS_Mesh* mesh = GetMeshDS();
   switch ( type ) {
-  case SMDSAbs_Edge:
-    if ( nbnode == 2 )
-      if ( ID ) e = mesh->AddEdgeWithID(node[0], node[1], ID);
-      else      e = mesh->AddEdge      (node[0], node[1] );
-    else if ( nbnode == 3 )
-      if ( ID ) e = mesh->AddEdgeWithID(node[0], node[1], node[2], ID);
-      else      e = mesh->AddEdge      (node[0], node[1], node[2] );
-    break;
   case SMDSAbs_Face:
     if ( !isPoly ) {
-      if      (nbnode == 3)
-        if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], ID);
-        else      e = mesh->AddFace      (node[0], node[1], node[2] );
-      else if (nbnode == 4) 
-        if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], ID);
-        else      e = mesh->AddFace      (node[0], node[1], node[2], node[3] );
-      else if (nbnode == 6)
-        if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
-                                          node[4], node[5], ID);
-        else      e = mesh->AddFace      (node[0], node[1], node[2], node[3],
-                                          node[4], node[5] );
-      else if (nbnode == 8)
-        if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
-                                          node[4], node[5], node[6], node[7], ID);
-        else      e = mesh->AddFace      (node[0], node[1], node[2], node[3],
-                                          node[4], node[5], node[6], node[7] );
+      if      (nbnode == 3) {
+        if ( ID >= 0 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], ID);
+        else           e = mesh->AddFace      (node[0], node[1], node[2] );
+      }
+      else if (nbnode == 4) {
+        if ( ID >= 0 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], ID);
+        else           e = mesh->AddFace      (node[0], node[1], node[2], node[3] );
+      }
+      else if (nbnode == 6) {
+        if ( ID >= 0 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
+                                              node[4], node[5], ID);
+        else           e = mesh->AddFace      (node[0], node[1], node[2], node[3],
+                                              node[4], node[5] );
+      }
+      else if (nbnode == 8) {
+        if ( ID >= 0 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
+                                              node[4], node[5], node[6], node[7], ID);
+        else           e = mesh->AddFace      (node[0], node[1], node[2], node[3],
+                                              node[4], node[5], node[6], node[7] );
+      }
     } else {
-      if ( ID ) e = mesh->AddPolygonalFaceWithID(node, ID);
-      else      e = mesh->AddPolygonalFace      (node    );
+      if ( ID >= 0 ) e = mesh->AddPolygonalFaceWithID(node, ID);
+      else           e = mesh->AddPolygonalFace      (node    );
     }
     break;
+
   case SMDSAbs_Volume:
     if ( !isPoly ) {
-      if      (nbnode == 4)
-        if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3] );
-      else if (nbnode == 5)
-        if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
-                                            node[4], ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
-                                            node[4] );
-      else if (nbnode == 6)
-        if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
-                                            node[4], node[5], ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
-                                            node[4], node[5] );
-      else if (nbnode == 8)
-        if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7], ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7] );
-      else if (nbnode == 10)
-        if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9], ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9] );
-      else if (nbnode == 13)
-        if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9], node[10],node[11],
-                                            node[12],ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9], node[10],node[11],
-                                            node[12] );
-      else if (nbnode == 15)
-        if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9], node[10],node[11],
-                                            node[12],node[13],node[14],ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9], node[10],node[11],
-                                            node[12],node[13],node[14] );
-      else if (nbnode == 20)
-        if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9], node[10],node[11],
-                                            node[12],node[13],node[14],node[15],
-                                            node[16],node[17],node[18],node[19],ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9], node[10],node[11],
-                                            node[12],node[13],node[14],node[15],
-                                            node[16],node[17],node[18],node[19] );
+      if      (nbnode == 4) {
+        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3] );
+      }
+      else if (nbnode == 5) {
+        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
+                                                node[4], ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
+                                                node[4] );
+      }
+      else if (nbnode == 6) {
+        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
+                                                node[4], node[5], ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
+                                                node[4], node[5] );
+      }
+      else if (nbnode == 8) {
+        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
+                                                node[4], node[5], node[6], node[7], ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
+                                                node[4], node[5], node[6], node[7] );
+      }
+      else if (nbnode == 10) {
+        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
+                                                node[4], node[5], node[6], node[7],
+                                                node[8], node[9], ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
+                                                node[4], node[5], node[6], node[7],
+                                                node[8], node[9] );
+      }
+      else if (nbnode == 13) {
+        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
+                                                node[4], node[5], node[6], node[7],
+                                                node[8], node[9], node[10],node[11],
+                                                node[12],ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
+                                                node[4], node[5], node[6], node[7],
+                                                node[8], node[9], node[10],node[11],
+                                                node[12] );
+      }
+      else if (nbnode == 15) {
+        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
+                                                node[4], node[5], node[6], node[7],
+                                                node[8], node[9], node[10],node[11],
+                                                node[12],node[13],node[14],ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
+                                                node[4], node[5], node[6], node[7],
+                                                node[8], node[9], node[10],node[11],
+                                                node[12],node[13],node[14] );
+      }
+      else if (nbnode == 20) {
+        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
+                                                node[4], node[5], node[6], node[7],
+                                                node[8], node[9], node[10],node[11],
+                                                node[12],node[13],node[14],node[15],
+                                                node[16],node[17],node[18],node[19],ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
+                                                node[4], node[5], node[6], node[7],
+                                                node[8], node[9], node[10],node[11],
+                                                node[12],node[13],node[14],node[15],
+                                                node[16],node[17],node[18],node[19] );
+      }
+    }
+    break;
+
+  case SMDSAbs_Edge:
+    if ( nbnode == 2 ) {
+      if ( ID >= 0 ) e = mesh->AddEdgeWithID(node[0], node[1], ID);
+      else           e = mesh->AddEdge      (node[0], node[1] );
+    }
+    else if ( nbnode == 3 ) {
+      if ( ID >= 0 ) e = mesh->AddEdgeWithID(node[0], node[1], node[2], ID);
+      else           e = mesh->AddEdge      (node[0], node[1], node[2] );
     }
+    break;
+
+  case SMDSAbs_0DElement:
+    if ( nbnode == 1 ) {
+      if ( ID >= 0 ) e = mesh->Add0DElementWithID(node[0], ID);
+      else           e = mesh->Add0DElement      (node[0] );
+    }
+    break;
+
+  case SMDSAbs_Node:
+    if ( ID >= 0 ) e = mesh->AddNodeWithID(node[0]->X(), node[0]->Y(), node[0]->Z(), ID);
+    else           e = mesh->AddNode      (node[0]->X(), node[0]->Y(), node[0]->Z());
+    break;
+
+  default:;
   }
+  if ( e ) myLastCreatedElems.Append( e );
   return e;
 }
 
@@ -258,8 +285,8 @@ SMDS_MeshElement* SMESH_MeshEditor::AddElement(const vector<int> &       nodeIDs
 //           Modify a compute state of sub-meshes which become empty
 //=======================================================================
 
-bool SMESH_MeshEditor::Remove (const list< int >& theIDs,
-                               const bool         isNodes )
+int SMESH_MeshEditor::Remove (const list< int >& theIDs,
+                              const bool         isNodes )
 {
   myLastCreatedElems.Clear();
   myLastCreatedNodes.Clear();
@@ -267,6 +294,7 @@ bool SMESH_MeshEditor::Remove (const list< int >& theIDs,
   SMESHDS_Mesh* aMesh = GetMeshDS();
   set< SMESH_subMesh *> smmap;
 
+  int removed = 0;
   list<int>::const_iterator it = theIDs.begin();
   for ( ; it != theIDs.end(); it++ ) {
     const SMDS_MeshElement * elem;
@@ -281,7 +309,7 @@ bool SMESH_MeshEditor::Remove (const list< int >& theIDs,
     if ( isNodes ) {
       const SMDS_MeshNode* node = cast2Node( elem );
       if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX )
-        if ( int aShapeID = node->GetPosition()->GetShapeId() )
+        if ( int aShapeID = node->getshapeId() )
           if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
             smmap.insert( sm );
     }
@@ -303,6 +331,7 @@ bool SMESH_MeshEditor::Remove (const list< int >& theIDs,
       aMesh->RemoveNode( static_cast< const SMDS_MeshNode* >( elem ));
     else
       aMesh->RemoveElement( elem );
+    removed++;
   }
 
   // Notify sub-meshes about modification
@@ -316,7 +345,7 @@ bool SMESH_MeshEditor::Remove (const list< int >& theIDs,
   //   if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( 1 ) )
   //     sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
 
-  return true;
+  return removed;
 }
 
 //=======================================================================
@@ -334,22 +363,21 @@ int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem)
   if ( aMesh->ShapeToMesh().IsNull() )
     return 0;
 
-  if ( theElem->GetType() == SMDSAbs_Node ) {
-    const SMDS_PositionPtr& aPosition =
-      static_cast<const SMDS_MeshNode*>( theElem )->GetPosition();
-    if ( aPosition.get() )
-      return aPosition->GetShapeId();
-    else
-      return 0;
-  }
+  if ( theElem->GetType() == SMDSAbs_Node )
+    {
+      int aShapeID = theElem->getshapeId();
+      if (aShapeID <= 0)
+        return 0;
+      else
+        return aShapeID;
+    }
 
   TopoDS_Shape aShape; // the shape a node is on
   SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
   while ( nodeIt->more() ) {
     const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
-    const SMDS_PositionPtr& aPosition = node->GetPosition();
-    if ( aPosition.get() ) {
-      int aShapeID = aPosition->GetShapeId();
+    int aShapeID = node->getshapeId();
+    if (aShapeID > 0) {
       SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID );
       if ( sm ) {
         if ( sm->Contains( theElem ))
@@ -480,15 +508,19 @@ static bool GetNodesFromTwoTria(const SMDS_MeshElement * theTria1,
 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1,
                                     const SMDS_MeshElement * theTria2 )
 {
+  MESSAGE("InverseDiag");
   myLastCreatedElems.Clear();
   myLastCreatedNodes.Clear();
 
   if (!theTria1 || !theTria2)
     return false;
 
-  const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( theTria1 );
-  const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( theTria2 );
-  if (F1 && F2) {
+  const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( theTria1 );
+  if (!F1) return false;
+  const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( theTria2 );
+  if (!F2) return false;
+  if ((theTria1->GetEntityType() == SMDSEntity_Triangle) &&
+      (theTria2->GetEntityType() == SMDSEntity_Triangle)) {
 
     //  1 +--+ A  theTria1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
     //    | /|    theTria2: ( B A 2 ) B->1 ( 1 A 2 )   |\ |
@@ -525,12 +557,14 @@ bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1,
     // find indices of 1,2 and of A,B in theTria1
     int iA = 0, iB = 0, i1 = 0, i2 = 0;
     for ( i = 0; i < 6; i++ ) {
-      if ( sameInd [ i ] == 0 )
+      if ( sameInd [ i ] == 0 ) {
         if ( i < 3 ) i1 = i;
         else         i2 = i;
-      else if (i < 3)
+      }
+      else if (i < 3) {
         if ( iA ) iB = i;
         else      iA = i;
+      }
     }
     // nodes 1 and 2 should not be the same
     if ( aNodes[ i1 ] == aNodes[ i2 ] )
@@ -541,24 +575,18 @@ bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1,
     // theTria2: B->1
     aNodes[ sameInd[ iB ]] = aNodes[ i1 ];
 
-    //MESSAGE( theTria1 << theTria2 );
-
     GetMeshDS()->ChangeElementNodes( theTria1, aNodes, 3 );
     GetMeshDS()->ChangeElementNodes( theTria2, &aNodes[ 3 ], 3 );
 
-    //MESSAGE( theTria1 << theTria2 );
-
     return true;
 
   } // end if(F1 && F2)
 
   // check case of quadratic faces
-  const SMDS_QuadraticFaceOfNodes* QF1 =
-    dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (theTria1);
-  if(!QF1) return false;
-  const SMDS_QuadraticFaceOfNodes* QF2 =
-    dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (theTria2);
-  if(!QF2) return false;
+  if (theTria1->GetEntityType() != SMDSEntity_Quad_Triangle)
+    return false;
+  if (theTria2->GetEntityType() != SMDSEntity_Quad_Triangle)
+    return false;
 
   //       5
   //  1 +--+--+ 2  theTria1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
@@ -623,7 +651,7 @@ static bool findTriangles(const SMDS_MeshNode *    theNode1,
   it = theNode2->GetInverseElementIterator(SMDSAbs_Face);
   while (it->more()) {
     const SMDS_MeshElement* elem = it->next();
-    if ( emap.find( elem ) != emap.end() )
+    if ( emap.find( elem ) != emap.end() ) {
       if ( theTria1 ) {
         // theTria1 must be element with minimum ID
         if( theTria1->GetID() < elem->GetID() ) {
@@ -638,6 +666,7 @@ static bool findTriangles(const SMDS_MeshNode *    theNode1,
       else {
         theTria1 = elem;
       }
+    }
   }
   return ( theTria1 && theTria2 );
 }
@@ -661,11 +690,12 @@ bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshNode * theNode1,
   if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
     return false;
 
-  const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( tr1 );
-  //if (!F1) return false;
-  const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( tr2 );
-  //if (!F2) return false;
-  if (F1 && F2) {
+  const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( tr1 );
+  if (!F1) return false;
+  const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( tr2 );
+  if (!F2) return false;
+  if ((tr1->GetEntityType() == SMDSEntity_Triangle) &&
+      (tr2->GetEntityType() == SMDSEntity_Triangle)) {
 
     //  1 +--+ A  tr1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
     //    | /|    tr2: ( B A 2 ) B->1 ( 1 A 2 )   |\ |
@@ -703,23 +733,13 @@ bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshNode * theNode1,
     // tr2: B->1
     aNodes2[ iB2 ] = aNodes1[ i1 ];
 
-    //MESSAGE( tr1 << tr2 );
-
     GetMeshDS()->ChangeElementNodes( tr1, aNodes1, 3 );
     GetMeshDS()->ChangeElementNodes( tr2, aNodes2, 3 );
 
-    //MESSAGE( tr1 << tr2 );
-
     return true;
   }
 
   // check case of quadratic faces
-  const SMDS_QuadraticFaceOfNodes* QF1 =
-    dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr1);
-  if(!QF1) return false;
-  const SMDS_QuadraticFaceOfNodes* QF2 =
-    dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr2);
-  if(!QF2) return false;
   return InverseDiag(tr1,tr2);
 }
 
@@ -793,34 +813,39 @@ bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
   if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
     return false;
 
-  const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( tr1 );
-  //if (!F1) return false;
-  const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( tr2 );
-  //if (!F2) return false;
-  if (F1 && F2) {
+  const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( tr1 );
+  if (!F1) return false;
+  const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( tr2 );
+  if (!F2) return false;
+  SMESHDS_Mesh * aMesh = GetMeshDS();
+
+  if ((tr1->GetEntityType() == SMDSEntity_Triangle) &&
+      (tr2->GetEntityType() == SMDSEntity_Triangle)) {
 
     const SMDS_MeshNode* aNodes [ 4 ];
     if ( ! getQuadrangleNodes( aNodes, theNode1, theNode2, tr1, tr2 ))
       return false;
 
-    //MESSAGE( endl << tr1 << tr2 );
-
-    GetMeshDS()->ChangeElementNodes( tr1, aNodes, 4 );
-    myLastCreatedElems.Append(tr1);
-    GetMeshDS()->RemoveElement( tr2 );
-
-    //MESSAGE( endl << tr1 );
+    const SMDS_MeshElement* newElem = 0;
+    newElem = aMesh->AddFace( aNodes[0], aNodes[1], aNodes[2], aNodes[3] );
+    myLastCreatedElems.Append(newElem);
+    AddToSameGroups( newElem, tr1, aMesh );
+    int aShapeId = tr1->getshapeId();
+    if ( aShapeId )
+      {
+        aMesh->SetMeshElementOnShape( newElem, aShapeId );
+      }
+    aMesh->RemoveElement( tr1 );
+    aMesh->RemoveElement( tr2 );
 
     return true;
   }
 
   // check case of quadratic faces
-  const SMDS_QuadraticFaceOfNodes* QF1 =
-    dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr1);
-  if(!QF1) return false;
-  const SMDS_QuadraticFaceOfNodes* QF2 =
-    dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr2);
-  if(!QF2) return false;
+  if (tr1->GetEntityType() != SMDSEntity_Quad_Triangle)
+    return false;
+  if (tr2->GetEntityType() != SMDSEntity_Quad_Triangle)
+    return false;
 
   //       5
   //  1 +--+--+ 2  tr1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
@@ -850,9 +875,18 @@ bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
   aNodes[6] = N2[3];
   aNodes[7] = N1[5];
 
-  GetMeshDS()->ChangeElementNodes( tr1, aNodes, 8 );
-  myLastCreatedElems.Append(tr1);
-  GetMeshDS()->RemoveElement( tr2 );
+  const SMDS_MeshElement* newElem = 0;
+  newElem = aMesh->AddFace( aNodes[0], aNodes[1], aNodes[2], aNodes[3],
+                            aNodes[4], aNodes[5], aNodes[6], aNodes[7]);
+  myLastCreatedElems.Append(newElem);
+  AddToSameGroups( newElem, tr1, aMesh );
+  int aShapeId = tr1->getshapeId();
+  if ( aShapeId )
+    {
+      aMesh->SetMeshElementOnShape( newElem, aShapeId );
+    }
+  aMesh->RemoveElement( tr1 );
+  aMesh->RemoveElement( tr2 );
 
   // remove middle node (9)
   GetMeshDS()->RemoveNode( N1[4] );
@@ -867,6 +901,7 @@ bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
 
 bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
 {
+  MESSAGE("Reorient");
   myLastCreatedElems.Clear();
   myLastCreatedNodes.Clear();
 
@@ -913,8 +948,10 @@ bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
   }
   case SMDSAbs_Volume: {
     if (theElem->IsPoly()) {
-      const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
-        static_cast<const SMDS_PolyhedralVolumeOfNodes*>( theElem );
+      // TODO reorient vtk polyhedron
+      MESSAGE("reorient vtk polyhedron ?");
+      const SMDS_VtkVolume* aPolyedre =
+        dynamic_cast<const SMDS_VtkVolume*>( theElem );
       if (!aPolyedre) {
         MESSAGE("Warning: bad volumic element");
         return false;
@@ -943,6 +980,7 @@ bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
       if ( !vTool.Set( theElem ))
         return false;
       vTool.Inverse();
+      MESSAGE("ChangeElementNodes reorient: check vTool.Inverse");
       return GetMeshDS()->ChangeElementNodes( theElem, vTool.GetNodes(), vTool.NbNodes() );
     }
   }
@@ -1015,21 +1053,21 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet &                   theElems,
     aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
 
     int aShapeId = FindShape( elem );
-    const SMDS_MeshElement* newElem = 0;
+    const SMDS_MeshElement* newElem1 = 0;
+    const SMDS_MeshElement* newElem2 = 0;
 
     if( !elem->IsQuadratic() ) {
 
       // split liner quadrangle
-
       if ( aBadRate1 <= aBadRate2 ) {
         // tr1 + tr2 is better
-        aMesh->ChangeElementNodes( elem, aNodes, 3 );
-        newElem = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
+        newElem1 = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
+        newElem2 = aMesh->AddFace( aNodes[2], aNodes[0], aNodes[1] );
       }
       else {
         // tr3 + tr4 is better
-        aMesh->ChangeElementNodes( elem, &aNodes[1], 3 );
-        newElem = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
+        newElem1 = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
+        newElem2 = aMesh->AddFace( aNodes[3], aNodes[1], aNodes[2] );
       }
     }
     else {
@@ -1090,8 +1128,10 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet &                   theElems,
         N[3] = aNodes[4];
         N[4] = aNodes[5];
         N[5] = newN;
-        newElem = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
-                                 aNodes[6], aNodes[7], newN );
+        newElem1 = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
+                                  aNodes[6], aNodes[7], newN );
+        newElem2 = aMesh->AddFace(aNodes[2], aNodes[0], aNodes[1],
+                                  newN,      aNodes[4], aNodes[5] );
       }
       else {
         N[0] = aNodes[1];
@@ -1100,21 +1140,27 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet &                   theElems,
         N[3] = aNodes[5];
         N[4] = aNodes[6];
         N[5] = newN;
-        newElem = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
-                                 aNodes[7], aNodes[4], newN );
+        newElem1 = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
+                                  aNodes[7], aNodes[4], newN );
+        newElem2 = aMesh->AddFace(aNodes[3], aNodes[1], aNodes[2],
+                                  newN,      aNodes[5], aNodes[6] );
       }
-      aMesh->ChangeElementNodes( elem, N, 6 );
-
     } // quadratic case
 
     // care of a new element
 
-    myLastCreatedElems.Append(newElem);
-    AddToSameGroups( newElem, elem, aMesh );
+    myLastCreatedElems.Append(newElem1);
+    myLastCreatedElems.Append(newElem2);
+    AddToSameGroups( newElem1, elem, aMesh );
+    AddToSameGroups( newElem2, elem, aMesh );
 
     // put a new triangle on the same shape
     if ( aShapeId )
-      aMesh->SetMeshElementOnShape( newElem, aShapeId );
+      {
+        aMesh->SetMeshElementOnShape( newElem1, aShapeId );
+        aMesh->SetMeshElementOnShape( newElem2, aShapeId );
+      }
+    aMesh->RemoveElement( elem );
   }
   return true;
 }
@@ -1123,6 +1169,7 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet &                   theElems,
 //function : BestSplit
 //purpose  : Find better diagonal for cutting.
 //=======================================================================
+
 int SMESH_MeshEditor::BestSplit (const SMDS_MeshElement*              theQuad,
                                  SMESH::Controls::NumericalFunctorPtr theCrit)
 {
@@ -1164,6 +1211,542 @@ int SMESH_MeshEditor::BestSplit (const SMDS_MeshElement*              theQuad,
   return -1;
 }
 
+namespace
+{
+  // Methods of splitting volumes into tetra
+
+  const int theHexTo5_1[5*4+1] =
+    {
+      0, 1, 2, 5,    0, 4, 5, 7,     0, 2, 3, 7,    2, 5, 6, 7,     0, 5, 2, 7,   -1
+    };
+  const int theHexTo5_2[5*4+1] =
+    {
+      1, 2, 3, 6,    1, 4, 5, 6,     0, 1, 3, 4,    3, 4, 6, 7,     1, 3, 4, 6,   -1
+    };
+  const int* theHexTo5[2] = { theHexTo5_1, theHexTo5_2 };
+
+  const int theHexTo6_1[6*4+1] =
+    {
+      1, 5, 6, 0,    0, 1, 2, 6,     0, 4, 5, 6,    0, 4, 6, 7,     0, 2, 3, 6,   0, 3, 7, 6,  -1
+    };
+  const int theHexTo6_2[6*4+1] =
+    {
+      2, 6, 7, 1,    1, 2, 3, 7,     1, 5, 6, 7,    1, 5, 7, 4,     1, 3, 0, 7,   1, 0, 4, 7,  -1
+    };
+  const int theHexTo6_3[6*4+1] =
+    {
+      3, 7, 4, 2,    2, 3, 0, 4,     2, 6, 7, 4,    2, 6, 4, 5,     2, 0, 1, 4,   2, 1, 5, 4,  -1
+    };
+  const int theHexTo6_4[6*4+1] =
+    {
+      0, 4, 5, 3,    3, 0, 1, 5,     3, 7, 4, 5,    3, 7, 5, 6,     3, 1, 2, 5,   3, 2, 6, 5,  -1
+    };
+  const int* theHexTo6[4] = { theHexTo6_1, theHexTo6_2, theHexTo6_3, theHexTo6_4 };
+
+  const int thePyraTo2_1[2*4+1] =
+    {
+      0, 1, 2, 4,    0, 2, 3, 4,   -1
+    };
+  const int thePyraTo2_2[2*4+1] =
+    {
+      1, 2, 3, 4,    1, 3, 0, 4,   -1
+    };
+  const int* thePyraTo2[2] = { thePyraTo2_1, thePyraTo2_2 };
+
+  const int thePentaTo3_1[3*4+1] =
+    {
+      0, 1, 2, 3,    1, 3, 4, 2,     2, 3, 4, 5,    -1
+    };
+  const int thePentaTo3_2[3*4+1] =
+    {
+      1, 2, 0, 4,    2, 4, 5, 0,     0, 4, 5, 3,    -1
+    };
+  const int thePentaTo3_3[3*4+1] =
+    {
+      2, 0, 1, 5,    0, 5, 3, 1,     1, 5, 3, 4,    -1
+    };
+  const int thePentaTo3_4[3*4+1] =
+    {
+      0, 1, 2, 3,    1, 3, 4, 5,     2, 3, 1, 5,    -1
+    };
+  const int thePentaTo3_5[3*4+1] =
+    {
+      1, 2, 0, 4,    2, 4, 5, 3,     0, 4, 2, 3,    -1
+    };
+  const int thePentaTo3_6[3*4+1] =
+    {
+      2, 0, 1, 5,    0, 5, 3, 4,     1, 5, 0, 4,    -1
+    };
+  const int* thePentaTo3[6] = { thePentaTo3_1, thePentaTo3_2, thePentaTo3_3,
+                                thePentaTo3_4, thePentaTo3_5, thePentaTo3_6 };
+
+  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;
+  };
+  struct TSplitMethod
+  {
+    int        _nbTetra;
+    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) {}
+    ~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;
+      return false;
+    }
+  };
+
+  //=======================================================================
+  /*!
+   * \brief return TSplitMethod for the given element
+   */
+  //=======================================================================
+
+  TSplitMethod getSplitMethod( SMDS_VolumeTool& vol, const int theMethodFlags)
+  {
+    const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
+
+    // at HEXA_TO_24 method, each face of volume is split into triangles each based on
+    // an edge and a face barycenter; tertaherdons are based on triangles and
+    // a volume barycenter
+    const bool is24TetMode = ( theMethodFlags == SMESH_MeshEditor::HEXA_TO_24 );
+
+    // Find out how adjacent volumes are split
+
+    vector < list< TTriangleFacet > > triaSplitsByFace( vol.NbFaces() ); // splits of each side
+    int hasAdjacentSplits = 0, maxTetConnSize = 0;
+    for ( int iF = 0; iF < vol.NbFaces(); ++iF )
+    {
+      int nbNodes = vol.NbFaceNodes( iF ) / iQ;
+      maxTetConnSize += 4 * ( nbNodes - (is24TetMode ? 0 : 2));
+      if ( nbNodes < 4 ) continue;
+
+      list< TTriangleFacet >& triaSplits = triaSplitsByFace[ iF ];
+      const int* nInd = vol.GetFaceNodesIndices( iF );
+      if ( nbNodes == 4 )
+      {
+        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 );
+      }
+      else
+      {
+        int iCom = 0; // common node of triangle faces to split into
+        for ( int iVar = 0; iVar < nbNodes; ++iVar, ++iCom )
+        {
+          TTriangleFacet t012( nInd[ iQ * ( iCom             )],
+                               nInd[ iQ * ( (iCom+1)%nbNodes )],
+                               nInd[ iQ * ( (iCom+2)%nbNodes )]);
+          TTriangleFacet t023( nInd[ iQ * ( iCom             )],
+                               nInd[ iQ * ( (iCom+2)%nbNodes )],
+                               nInd[ iQ * ( (iCom+3)%nbNodes )]);
+          if ( t012.hasAdjacentTetra( vol.Element() ) && t023.hasAdjacentTetra( vol.Element() ))
+          {
+            triaSplits.push_back( t012 );
+            triaSplits.push_back( t023 );
+            break;
+          }
+        }
+      }
+      if ( !triaSplits.empty() )
+        hasAdjacentSplits = true;
+    }
+
+    // Among variants of split method select one compliant with adjacent volumes
+
+    TSplitMethod method;
+    if ( !vol.Element()->IsPoly() && !is24TetMode )
+    {
+      int nbVariants = 2, nbTet = 0;
+      const int** connVariants = 0;
+      switch ( vol.Element()->GetEntityType() )
+      {
+      case SMDSEntity_Hexa:
+      case SMDSEntity_Quad_Hexa:
+        if ( theMethodFlags == SMESH_MeshEditor::HEXA_TO_5 )
+          connVariants = theHexTo5, nbTet = 5;
+        else
+          connVariants = theHexTo6, nbTet = 6, nbVariants = 4;
+        break;
+      case SMDSEntity_Pyramid:
+      case SMDSEntity_Quad_Pyramid:
+        connVariants = thePyraTo2;  nbTet = 2;
+        break;
+      case SMDSEntity_Penta:
+      case SMDSEntity_Quad_Penta:
+        connVariants = thePentaTo3; nbTet = 3; nbVariants = 6;
+        break;
+      default:
+        nbVariants = 0;
+      }
+      for ( int variant = 0; variant < nbVariants && method._nbTetra == 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 )
+        {
+          bool facetCreated = true;
+          for ( int iF = 0; facetCreated && iF < triaSplitsByFace.size(); ++iF )
+          {
+            list< TTriangleFacet >::const_iterator facet = triaSplitsByFace[iF].begin();
+            for ( ; facetCreated && facet != triaSplitsByFace[iF].end(); ++facet )
+              facetCreated = method.hasFacet( *facet );
+          }
+          if ( !facetCreated )
+            method = TSplitMethod(0); // incompatible method
+        }
+      }
+    }
+    if ( method._nbTetra < 1 )
+    {
+      // No standard method is applicable, use a generic solution:
+      // each facet of a volume is split into triangles and
+      // each of triangles and a volume barycenter form a tetrahedron.
+
+      int* connectivity = new int[ maxTetConnSize + 1 ];
+      method._connectivity = connectivity;
+      method._ownConn = true;
+      method._baryNode = true;
+
+      int connSize = 0;
+      int baryCenInd = vol.NbNodes();
+      for ( int iF = 0; iF < vol.NbFaces(); ++iF )
+      {
+        const int nbNodes = vol.NbFaceNodes( iF ) / iQ;
+        const int*   nInd = vol.GetFaceNodesIndices( iF );
+        // find common node of triangle facets of tetra to create
+        int iCommon = 0; // index in linear numeration
+        const list< TTriangleFacet >& triaSplits = triaSplitsByFace[ iF ];
+        if ( !triaSplits.empty() )
+        {
+          // by found facets
+          const TTriangleFacet* facet = &triaSplits.front();
+          for ( ; iCommon < nbNodes-1 ; ++iCommon )
+            if ( facet->contains( nInd[ iQ * iCommon ]) &&
+                 facet->contains( nInd[ iQ * ((iCommon+2)%nbNodes) ]))
+              break;
+        }
+        else if ( nbNodes > 3 && !is24TetMode )
+        {
+          // find the best method of splitting into triangles by aspect ratio
+          SMESH::Controls::NumericalFunctorPtr aspectRatio( new SMESH::Controls::AspectRatio);
+          map< double, int > badness2iCommon;
+          const SMDS_MeshNode** nodes = vol.GetFaceNodes( iF );
+          int nbVariants = ( nbNodes == 4 ? 2 : nbNodes );
+          for ( int iVar = 0; iVar < nbVariants; ++iVar, ++iCommon )
+            for ( int iLast = iCommon+2; iLast < iCommon+nbNodes; ++iLast )
+            {
+              SMDS_FaceOfNodes tria ( nodes[ iQ*( iCommon         )],
+                                      nodes[ iQ*((iLast-1)%nbNodes)],
+                                      nodes[ iQ*((iLast  )%nbNodes)]);
+              double badness = getBadRate( &tria, aspectRatio );
+              badness2iCommon.insert( make_pair( badness, iCommon ));
+            }
+          // use iCommon with lowest badness
+          iCommon = badness2iCommon.begin()->second;
+        }
+        if ( iCommon >= nbNodes )
+          iCommon = 0; // something wrong
+
+        // fill connectivity of tetrahedra based on a current face
+        int nbTet = nbNodes - 2;
+        if ( is24TetMode && nbNodes > 3 && triaSplits.empty())
+        {
+          method._faceBaryNode.insert( make_pair( iF, (const SMDS_MeshNode*)0 ));
+          int faceBaryCenInd = baryCenInd + method._faceBaryNode.size();
+          nbTet = nbNodes;
+          for ( int i = 0; i < nbTet; ++i )
+          {
+            int i1 = i, i2 = (i+1) % nbNodes;
+            if ( !vol.IsFaceExternal( iF )) swap( i1, i2 );
+            connectivity[ connSize++ ] = nInd[ iQ * i1 ];
+            connectivity[ connSize++ ] = nInd[ iQ * i2 ];
+            connectivity[ connSize++ ] = faceBaryCenInd;
+            connectivity[ connSize++ ] = baryCenInd;
+          }
+        }
+        else
+        {
+          for ( int i = 0; i < nbTet; ++i )
+          {
+            int i1 = (iCommon+1+i) % nbNodes, i2 = (iCommon+2+i) % nbNodes;
+            if ( !vol.IsFaceExternal( iF )) swap( i1, i2 );
+            connectivity[ connSize++ ] = nInd[ iQ * iCommon ];
+            connectivity[ connSize++ ] = nInd[ iQ * i1 ];
+            connectivity[ connSize++ ] = nInd[ iQ * i2 ];
+            connectivity[ connSize++ ] = baryCenInd;
+          }
+        }
+        method._nbTetra += nbTet;
+      }
+      connectivity[ connSize++ ] = -1;
+    }
+    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
+  {
+    // find the tetrahedron including the three nodes of facet
+    const SMDS_MeshNode* n1 = elem->GetNode(_n1);
+    const SMDS_MeshNode* n2 = elem->GetNode(_n2);
+    const SMDS_MeshNode* n3 = elem->GetNode(_n3);
+    SMDS_ElemIteratorPtr volIt1 = n1->GetInverseElementIterator(SMDSAbs_Volume);
+    while ( volIt1->more() )
+    {
+      const SMDS_MeshElement* v = volIt1->next();
+      if ( v->GetEntityType() != ( v->IsQuadratic() ? SMDSEntity_Quad_Tetra : SMDSEntity_Tetra ))
+        continue;
+      SMDS_ElemIteratorPtr volIt2 = n2->GetInverseElementIterator(SMDSAbs_Volume);
+      while ( volIt2->more() )
+        if ( v != volIt2->next() )
+          continue;
+      SMDS_ElemIteratorPtr volIt3 = n3->GetInverseElementIterator(SMDSAbs_Volume);
+      while ( volIt3->more() )
+        if ( v == volIt3->next() )
+          return true;
+    }
+    return false;
+  }
+
+  //=======================================================================
+  /*!
+   * \brief A key of a face of volume
+   */
+  //=======================================================================
+
+  struct TVolumeFaceKey: pair< int, pair< int, int> >
+  {
+    TVolumeFaceKey( SMDS_VolumeTool& vol, int iF )
+    {
+      TIDSortedNodeSet sortedNodes;
+      const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
+      int nbNodes = vol.NbFaceNodes( iF );
+      const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iF );
+      for ( int i = 0; i < nbNodes; i += iQ )
+        sortedNodes.insert( fNodes[i] );
+      TIDSortedNodeSet::iterator n = sortedNodes.begin();
+      first = (*(n++))->GetID();
+      second.first = (*(n++))->GetID();
+      second.second = (*(n++))->GetID();
+    }
+  };
+} // namespace
+
+//=======================================================================
+//function : SplitVolumesIntoTetra
+//purpose  : Split volumic elements into tetrahedra.
+//=======================================================================
+
+void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
+                                              const int                theMethodFlags)
+{
+  // std-like iterator on coordinates of nodes of mesh element
+  typedef SMDS_StdIterator< TNodeXYZ, SMDS_ElemIteratorPtr > NXyzIterator;
+  NXyzIterator xyzEnd;
+
+  SMDS_VolumeTool    volTool;
+  SMESH_MesherHelper helper( *GetMesh());
+
+  SMESHDS_SubMesh* subMesh = GetMeshDS()->MeshElements(1);
+  SMESHDS_SubMesh* fSubMesh = subMesh;
+  
+  SMESH_SequenceOfElemPtr newNodes, newElems;
+
+  // map face of volume to it's baricenrtic node
+  map< TVolumeFaceKey, const SMDS_MeshNode* > volFace2BaryNode;
+  double bc[3];
+
+  TIDSortedElemSet::const_iterator elem = theElems.begin();
+  for ( ; elem != theElems.end(); ++elem )
+  {
+    SMDSAbs_EntityType geomType = (*elem)->GetEntityType();
+    if ( geomType <= SMDSEntity_Quad_Tetra )
+      continue; // tetra or face or ...
+
+    if ( !volTool.Set( *elem )) continue; // not volume? strange...
+
+    TSplitMethod splitMethod = getSplitMethod( volTool, theMethodFlags );
+    if ( splitMethod._nbTetra < 1 ) continue;
+
+    // find submesh to add new tetras to
+    if ( !subMesh || !subMesh->Contains( *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() )
+    {
+      iQ = 2;
+      // add quadratic links to the helper
+      for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
+      {
+        const SMDS_MeshNode** fNodes = volTool.GetFaceNodes( iF );
+        for ( int iN = 0; iN < volTool.NbFaceNodes( iF ); iN += iQ )
+          helper.AddTLinkNode( fNodes[iF], fNodes[iF+2], fNodes[iF+1] );
+      }
+      helper.SetIsQuadratic( true );
+    }
+    else
+    {
+      iQ = 1;
+      helper.SetIsQuadratic( false );
+    }
+    vector<const SMDS_MeshNode*> nodes( (*elem)->begin_nodes(), (*elem)->end_nodes() );
+    if ( splitMethod._baryNode )
+    {
+      // make a node at barycenter
+      volTool.GetBaryCenter( bc[0], bc[1], bc[2] );
+      SMDS_MeshNode* gcNode = helper.AddNode( bc[0], bc[1], bc[2] );
+      nodes.push_back( gcNode );
+      newNodes.Append( gcNode );
+    }
+    if ( !splitMethod._faceBaryNode.empty() )
+    {
+      // make or find baricentric nodes of faces
+      map<int, const SMDS_MeshNode*>::iterator iF_n = splitMethod._faceBaryNode.begin();
+      for ( ; iF_n != splitMethod._faceBaryNode.end(); ++iF_n )
+      {
+        map< TVolumeFaceKey, const SMDS_MeshNode* >::iterator f_n =
+          volFace2BaryNode.insert
+          ( make_pair( TVolumeFaceKey( volTool,iF_n->first ), (const SMDS_MeshNode*)0) ).first;
+        if ( !f_n->second )
+        {
+          volTool.GetFaceBaryCenter( iF_n->first, bc[0], bc[1], bc[2] );
+          newNodes.Append( f_n->second = helper.AddNode( bc[0], bc[1], bc[2] ));
+        }
+        nodes.push_back( iF_n->second = f_n->second );
+      }
+    }
+
+    // make tetras
+    helper.SetElementsOnShape( true );
+    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] ]));
+
+    ReplaceElemInGroups( *elem, tetras, GetMeshDS() );
+
+    // Split faces on sides of the split volume
+
+    const SMDS_MeshNode** volNodes = volTool.GetNodes();
+    for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
+    {
+      const int nbNodes = volTool.NbFaceNodes( iF ) / iQ;
+      if ( nbNodes < 4 ) continue;
+
+      // find an existing face
+      vector<const SMDS_MeshNode*> fNodes( volTool.GetFaceNodes( iF ),
+                                           volTool.GetFaceNodes( iF ) + nbNodes*iQ );
+      while ( const SMDS_MeshElement* face = GetMeshDS()->FindFace( fNodes ))
+      {
+        // make triangles
+        helper.SetElementsOnShape( false );
+        vector< const SMDS_MeshElement* > triangles;
+
+        map<int, const SMDS_MeshNode*>::iterator iF_n = splitMethod._faceBaryNode.find(iF);
+        if ( iF_n != splitMethod._faceBaryNode.end() )
+        {
+          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;
+            if ( !volTool.IsFaceExternal( iF ))
+              swap( n2, n3 );
+            triangles.push_back( helper.AddFace( n1,n2,n3 ));
+          }
+        }
+        else
+        {
+          // among possible triangles create ones discribed by split method
+          const int* nInd = volTool.GetFaceNodesIndices( iF );
+          int nbVariants = ( nbNodes == 4 ? 2 : nbNodes );
+          int iCom = 0; // common node of triangle faces to split into
+          list< TTriangleFacet > facets;
+          for ( int iVar = 0; iVar < nbVariants; ++iVar, ++iCom )
+          {
+            TTriangleFacet t012( nInd[ iQ * ( iCom                )],
+                                 nInd[ iQ * ( (iCom+1)%nbNodes )],
+                                 nInd[ iQ * ( (iCom+2)%nbNodes )]);
+            TTriangleFacet t023( nInd[ iQ * ( iCom                )],
+                                 nInd[ iQ * ( (iCom+2)%nbNodes )],
+                                 nInd[ iQ * ( (iCom+3)%nbNodes )]);
+            if ( splitMethod.hasFacet( t012 ) && splitMethod.hasFacet( t023 ))
+            {
+              facets.push_back( t012 );
+              facets.push_back( t023 );
+              for ( int iLast = iCom+4; iLast < iCom+nbNodes; ++iLast )
+                facets.push_back( TTriangleFacet( nInd[ iQ * ( iCom             )],
+                                                  nInd[ iQ * ((iLast-1)%nbNodes )],
+                                                  nInd[ iQ * ((iLast  )%nbNodes )]));
+              break;
+            }
+          }
+          list< TTriangleFacet >::iterator facet = facets.begin();
+          for ( ; facet != facets.end(); ++facet )
+          {
+            if ( !volTool.IsFaceExternal( iF ))
+              swap( facet->_n2, facet->_n3 );
+            triangles.push_back( helper.AddFace( volNodes[ facet->_n1 ],
+                                                 volNodes[ facet->_n2 ],
+                                                 volNodes[ facet->_n3 ]));
+          }
+        }
+        // find submesh to add new triangles in
+        if ( !fSubMesh || !fSubMesh->Contains( face ))
+        {
+          int shapeID = FindShape( face );
+          fSubMesh = GetMeshDS()->MeshElements( shapeID );
+        }
+        for ( int i = 0; i < triangles.size(); ++i )
+        {
+          if ( !triangles[i] ) continue;
+          if ( fSubMesh )
+            fSubMesh->AddElement( triangles[i]);
+          newElems.Append( triangles[i] );
+        }
+        ReplaceElemInGroups( face, triangles, GetMeshDS() );
+        GetMeshDS()->RemoveFreeElement( face, fSubMesh, /*fromGroups=*/false );
+      }
+
+    } // loop on volume faces to split them into triangles
+
+    GetMeshDS()->RemoveFreeElement( *elem, subMesh, /*fromGroups=*/false );
+
+  } // loop on volumes to split
+
+  myLastCreatedNodes = newNodes;
+  myLastCreatedElems = newElems;
+}
+
 //=======================================================================
 //function : AddToSameGroups
 //purpose  : add elemToAdd to the groups the elemInGroups belongs to
@@ -1205,10 +1788,11 @@ void SMESH_MeshEditor::RemoveElemFromGroups (const SMDS_MeshElement* removeelem,
   }
 }
 
-//=======================================================================
-//function : ReplaceElemInGroups
-//purpose  : replace elemToRm by elemToAdd in the all groups
-//=======================================================================
+//================================================================================
+/*!
+ * \brief Replace elemToRm by elemToAdd in the all groups
+ */
+//================================================================================
 
 void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement* elemToRm,
                                             const SMDS_MeshElement* elemToAdd,
@@ -1225,6 +1809,29 @@ void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement* elemToRm,
   }
 }
 
+//================================================================================
+/*!
+ * \brief Replace elemToRm by elemToAdd in the all groups
+ */
+//================================================================================
+
+void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement*                elemToRm,
+                                            const vector<const SMDS_MeshElement*>& elemToAdd,
+                                            SMESHDS_Mesh *                         aMesh)
+{
+  const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
+  if (!groups.empty())
+  {
+    set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
+    for ( ; grIt != groups.end(); grIt++ ) {
+      SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
+      if ( group && group->SMDSGroup().Remove( elemToRm ) )
+        for ( int i = 0; i < elemToAdd.size(); ++i )
+          group->SMDSGroup().Add( elemToAdd[ i ] );
+    }
+  }
+}
+
 //=======================================================================
 //function : QuadToTri
 //purpose  : Cut quadrangles into triangles.
@@ -1261,20 +1868,28 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
         aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
 
       int aShapeId = FindShape( elem );
-      const SMDS_MeshElement* newElem = 0;
+      const SMDS_MeshElement* newElem1 = 0;
+      const SMDS_MeshElement* newElem2 = 0;
       if ( the13Diag ) {
-        aMesh->ChangeElementNodes( elem, aNodes, 3 );
-        newElem = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
+        newElem1 = aMesh->AddFace( aNodes[2], aNodes[0], aNodes[1] );
+        newElem2 = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
       }
       else {
-        aMesh->ChangeElementNodes( elem, &aNodes[1], 3 );
-        newElem = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
+        newElem1 = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
+        newElem2 = aMesh->AddFace( aNodes[3], aNodes[1], aNodes[2] );
       }
-      myLastCreatedElems.Append(newElem);
+      myLastCreatedElems.Append(newElem1);
+      myLastCreatedElems.Append(newElem2);
       // put a new triangle on the same shape and add to the same groups
       if ( aShapeId )
-        aMesh->SetMeshElementOnShape( newElem, aShapeId );
-      AddToSameGroups( newElem, elem, aMesh );
+        {
+          aMesh->SetMeshElementOnShape( newElem1, aShapeId );
+          aMesh->SetMeshElementOnShape( newElem2, aShapeId );
+        }
+      AddToSameGroups( newElem1, elem, aMesh );
+      AddToSameGroups( newElem2, elem, aMesh );
+      //aMesh->RemoveFreeElement(elem, aMesh->MeshElements(aShapeId), true);
+      aMesh->RemoveElement( elem );
     }
 
     // Quadratic quadrangle
@@ -1329,7 +1944,8 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
       myLastCreatedNodes.Append(newN);
 
       // create a new element
-      const SMDS_MeshElement* newElem = 0;
+      const SMDS_MeshElement* newElem1 = 0;
+      const SMDS_MeshElement* newElem2 = 0;
       const SMDS_MeshNode* N[6];
       if ( the13Diag ) {
         N[0] = aNodes[0];
@@ -1338,8 +1954,10 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
         N[3] = aNodes[4];
         N[4] = aNodes[5];
         N[5] = newN;
-        newElem = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
-                                 aNodes[6], aNodes[7], newN );
+        newElem1 = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
+                                  aNodes[6], aNodes[7], newN );
+        newElem2 = aMesh->AddFace(aNodes[2], aNodes[0], aNodes[1],
+                                  newN,      aNodes[4], aNodes[5] );
       }
       else {
         N[0] = aNodes[1];
@@ -1348,15 +1966,22 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
         N[3] = aNodes[5];
         N[4] = aNodes[6];
         N[5] = newN;
-        newElem = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
-                                 aNodes[7], aNodes[4], newN );
+        newElem1 = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
+                                  aNodes[7], aNodes[4], newN );
+        newElem2 = aMesh->AddFace(aNodes[3], aNodes[1], aNodes[2],
+                                  newN,      aNodes[5], aNodes[6] );
       }
-      myLastCreatedElems.Append(newElem);
-      aMesh->ChangeElementNodes( elem, N, 6 );
+      myLastCreatedElems.Append(newElem1);
+      myLastCreatedElems.Append(newElem2);
       // put a new triangle on the same shape and add to the same groups
       if ( aShapeId )
-        aMesh->SetMeshElementOnShape( newElem, aShapeId );
-      AddToSameGroups( newElem, elem, aMesh );
+        {
+          aMesh->SetMeshElementOnShape( newElem1, aShapeId );
+          aMesh->SetMeshElementOnShape( newElem2, aShapeId );
+        }
+      AddToSameGroups( newElem1, elem, aMesh );
+      AddToSameGroups( newElem2, elem, aMesh );
+      aMesh->RemoveElement( elem );
     }
   }
 
@@ -1402,7 +2027,7 @@ double getAngle(const SMDS_MeshElement * tr1,
     int i = 0, iDiag = -1;
     while ( it->more()) {
       const SMDS_MeshElement *n = it->next();
-      if ( n == n1 || n == n2 )
+      if ( n == n1 || n == n2 ) {
         if ( iDiag < 0)
           iDiag = i;
         else {
@@ -1412,6 +2037,7 @@ double getAngle(const SMDS_MeshElement * tr1,
             nFirst[ t ] = n;
           break;
         }
+      }
       i++;
     }
   }
@@ -1652,16 +2278,17 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
           mapEl_setLi.erase( tr2 );
           mapLi_listEl.erase( *link12 );
           if(tr1->NbNodes()==3) {
-            if( tr1->GetID() < tr2->GetID() ) {
-              aMesh->ChangeElementNodes( tr1, n12, 4 );
-              myLastCreatedElems.Append(tr1);
-              aMesh->RemoveElement( tr2 );
-            }
-            else {
-              aMesh->ChangeElementNodes( tr2, n12, 4 );
-              myLastCreatedElems.Append(tr2);
-              aMesh->RemoveElement( tr1);
-            }
+            const SMDS_MeshElement* newElem = 0;
+            newElem = aMesh->AddFace(n12[0], n12[1], n12[2], n12[3] );
+            myLastCreatedElems.Append(newElem);
+            AddToSameGroups( newElem, tr1, aMesh );
+            int aShapeId = tr1->getshapeId();
+            if ( aShapeId )
+              {
+                aMesh->SetMeshElementOnShape( newElem, aShapeId );
+              }
+            aMesh->RemoveElement( tr1 );
+            aMesh->RemoveElement( tr2 );
           }
           else {
             const SMDS_MeshNode* N1 [6];
@@ -1679,16 +2306,18 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
             aNodes[5] = N2[5];
             aNodes[6] = N2[3];
             aNodes[7] = N1[5];
-            if( tr1->GetID() < tr2->GetID() ) {
-              GetMeshDS()->ChangeElementNodes( tr1, aNodes, 8 );
-              myLastCreatedElems.Append(tr1);
-              GetMeshDS()->RemoveElement( tr2 );
-            }
-            else {
-              GetMeshDS()->ChangeElementNodes( tr2, aNodes, 8 );
-              myLastCreatedElems.Append(tr2);
-              GetMeshDS()->RemoveElement( tr1 );
-            }
+            const SMDS_MeshElement* newElem = 0;
+            newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3],
+                                     aNodes[4], aNodes[5], aNodes[6], aNodes[7]);
+            myLastCreatedElems.Append(newElem);
+            AddToSameGroups( newElem, tr1, aMesh );
+            int aShapeId = tr1->getshapeId();
+            if ( aShapeId )
+              {
+                aMesh->SetMeshElementOnShape( newElem, aShapeId );
+              }
+            aMesh->RemoveElement( tr1 );
+            aMesh->RemoveElement( tr2 );
             // remove middle node (9)
             GetMeshDS()->RemoveNode( N1[4] );
           }
@@ -1697,16 +2326,17 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
           mapEl_setLi.erase( tr3 );
           mapLi_listEl.erase( *link13 );
           if(tr1->NbNodes()==3) {
-            if( tr1->GetID() < tr2->GetID() ) {
-              aMesh->ChangeElementNodes( tr1, n13, 4 );
-              myLastCreatedElems.Append(tr1);
-              aMesh->RemoveElement( tr3 );
-            }
-            else {
-              aMesh->ChangeElementNodes( tr3, n13, 4 );
-              myLastCreatedElems.Append(tr3);
-              aMesh->RemoveElement( tr1 );
-            }
+            const SMDS_MeshElement* newElem = 0;
+            newElem = aMesh->AddFace(n13[0], n13[1], n13[2], n13[3] );
+            myLastCreatedElems.Append(newElem);
+            AddToSameGroups( newElem, tr1, aMesh );
+            int aShapeId = tr1->getshapeId();
+            if ( aShapeId )
+              {
+                aMesh->SetMeshElementOnShape( newElem, aShapeId );
+              }
+            aMesh->RemoveElement( tr1 );
+            aMesh->RemoveElement( tr3 );
           }
           else {
             const SMDS_MeshNode* N1 [6];
@@ -1724,16 +2354,18 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
             aNodes[5] = N2[5];
             aNodes[6] = N2[3];
             aNodes[7] = N1[5];
-            if( tr1->GetID() < tr2->GetID() ) {
-              GetMeshDS()->ChangeElementNodes( tr1, aNodes, 8 );
-              myLastCreatedElems.Append(tr1);
-              GetMeshDS()->RemoveElement( tr3 );
-            }
-            else {
-              GetMeshDS()->ChangeElementNodes( tr3, aNodes, 8 );
-              myLastCreatedElems.Append(tr3);
-              GetMeshDS()->RemoveElement( tr1 );
-            }
+            const SMDS_MeshElement* newElem = 0;
+            newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3],
+                                     aNodes[4], aNodes[5], aNodes[6], aNodes[7]);
+            myLastCreatedElems.Append(newElem);
+            AddToSameGroups( newElem, tr1, aMesh );
+            int aShapeId = tr1->getshapeId();
+            if ( aShapeId )
+              {
+                aMesh->SetMeshElementOnShape( newElem, aShapeId );
+              }
+            aMesh->RemoveElement( tr1 );
+            aMesh->RemoveElement( tr3 );
             // remove middle node (9)
             GetMeshDS()->RemoveNode( N1[4] );
           }
@@ -2345,7 +2977,7 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
       while ( nn++ < nbn ) {
         node = static_cast<const SMDS_MeshNode*>( itN->next() );
         const SMDS_PositionPtr& pos = node->GetPosition();
-        posType = pos.get() ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
+        posType = pos ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
         if (posType != SMDS_TOP_EDGE &&
             posType != SMDS_TOP_VERTEX &&
             theFixedNodes.find( node ) == theFixedNodes.end())
@@ -2403,27 +3035,27 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
         node = *n;
         gp_XY uv( 0, 0 );
         const SMDS_PositionPtr& pos = node->GetPosition();
-        posType = pos.get() ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
+        posType = pos ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
         // get existing UV
         switch ( posType ) {
         case SMDS_TOP_FACE: {
-          SMDS_FacePosition* fPos = ( SMDS_FacePosition* ) pos.get();
+          SMDS_FacePosition* fPos = ( SMDS_FacePosition* ) pos;
           uv.SetCoord( fPos->GetUParameter(), fPos->GetVParameter() );
           break;
         }
         case SMDS_TOP_EDGE: {
-          TopoDS_Shape S = aMesh->IndexToShape( pos->GetShapeId() );
+          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.get() )->GetUParameter();
+            double u = (( SMDS_EdgePosition* ) pos )->GetUParameter();
             uv = pcurve->Value( u ).XY();
           }
           break;
         }
         case SMDS_TOP_VERTEX: {
-          TopoDS_Shape S = aMesh->IndexToShape( pos->GetShapeId() );
+          TopoDS_Shape S = aMesh->IndexToShape( node->getshapeId() );
           if ( !S.IsNull() && S.ShapeType() == TopAbs_VERTEX )
             uv = BRep_Tool::Parameters( TopoDS::Vertex( S ), face ).XY();
           break;
@@ -2686,7 +3318,7 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
       if ( node_uv != uvMap.end() ) {
         gp_XY* uv = node_uv->second;
         node->SetPosition
-          ( SMDS_PositionPtr( new SMDS_FacePosition( *fId, uv->X(), uv->Y() )));
+          ( SMDS_PositionPtr( new SMDS_FacePosition( uv->X(), uv->Y() )));
       }
     }
 
@@ -2698,14 +3330,14 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
         helper.SetSubShape( face );
       list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
       for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
-        const SMDS_QuadraticFaceOfNodes* QF =
-          dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (*elemIt);
-        if(QF) {
+        const SMDS_VtkFace* QF =
+          dynamic_cast<const SMDS_VtkFace*> (*elemIt);
+        if(QF && QF->IsQuadratic()) {
           vector<const SMDS_MeshNode*> Ns;
           Ns.reserve(QF->NbNodes()+1);
-          SMDS_NodeIteratorPtr anIter = QF->interlacedNodesIterator();
+          SMDS_ElemIteratorPtr anIter = QF->interlacedNodesElemIterator();
           while ( anIter->more() )
-            Ns.push_back( anIter->next() );
+            Ns.push_back( cast2Node(anIter->next()) );
           Ns.push_back( Ns[0] );
           double x, y, z;
           for(int i=0; i<QF->NbNodes(); i=i+2) {
@@ -2783,6 +3415,7 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
                                     const int                             nbSteps,
                                     SMESH_SequenceOfElemPtr&              srcElements)
 {
+  //MESSAGE("sweepElement " << nbSteps);
   SMESHDS_Mesh* aMesh = GetMeshDS();
 
   // Loop on elem nodes:
@@ -2805,12 +3438,12 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
       return;
     }
 
-    issimple[iNode] = (listNewNodes.size()==nbSteps);
+    issimple[iNode] = (listNewNodes.size()==nbSteps); // is node medium
 
     itNN[ iNode ] = listNewNodes.begin();
     prevNod[ iNode ] = node;
     nextNod[ iNode ] = listNewNodes.front();
-    if( !issimple[iNode] ) {
+    if( !elem->IsQuadratic() || !issimple[iNode] ) {
       if ( prevNod[ iNode ] != nextNod [ iNode ])
         iNotSameNode = iNode;
       else {
@@ -2821,10 +3454,10 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
     }
   }
 
-  //cout<<"  nbSame = "<<nbSame<<endl;
+  //cerr<<"  nbSame = "<<nbSame<<endl;
   if ( nbSame == nbNodes || nbSame > 2) {
-    //MESSAGE( " Too many same nodes of element " << elem->GetID() );
-    INFOS( " Too many same nodes of element " << elem->GetID() );
+    MESSAGE( " Too many same nodes of element " << elem->GetID() );
+    //INFOS( " Too many same nodes of element " << elem->GetID() );
     return;
   }
 
@@ -2858,6 +3491,7 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
   }
 
   // make new elements
+  const SMDS_MeshElement* lastElem = elem;
   for (int iStep = 0; iStep < nbSteps; iStep++ ) {
     // get next nodes
     for ( iNode = 0; iNode < nbNodes; iNode++ ) {
@@ -2873,7 +3507,7 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
           nextNod[ iNode ] = *itNN[ iNode ];
           itNN[ iNode ]++;
         }
-        else if(!elem->IsQuadratic() || elem->IsMediumNode(prevNod[iNode]) ) {
+        else if(!elem->IsQuadratic() || lastElem->IsMediumNode(prevNod[iNode]) ) {
           // we have to use each second node
           //itNN[ iNode ]++;
           nextNod[ iNode ] = *itNN[ iNode ];
@@ -3178,6 +3812,7 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement*               elem,
       newElems.push_back( aNewElem );
       myLastCreatedElems.Append(aNewElem);
       srcElements.Append( elem );
+      lastElem = aNewElem;
     }
 
     // set new prev nodes
@@ -3206,6 +3841,7 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
                                   const int                nbSteps,
                                   SMESH_SequenceOfElemPtr& srcElements)
 {
+  MESSAGE("makeWalls");
   ASSERT( newElemsMap.size() == elemNewNodesMap.size() );
   SMESHDS_Mesh* aMesh = GetMeshDS();
 
@@ -3406,7 +4042,10 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
               if ( !f )
                 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ));
               else if ( nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
-                aMesh->ChangeElementNodes( f, nodes, nbn );
+                {
+                  myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ));
+                  aMesh->RemoveElement(f);
+                }
               break;
             }
             case 4: { ///// quadrangle
@@ -3414,7 +4053,10 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
               if ( !f )
                 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ));
               else if ( nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
-                aMesh->ChangeElementNodes( f, nodes, nbn );
+                {
+                  myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ));
+                  aMesh->RemoveElement(f);
+                }
               break;
             }
             default:
@@ -3434,7 +4076,9 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
                     tmpnodes[3] = nodes[1];
                     tmpnodes[4] = nodes[3];
                     tmpnodes[5] = nodes[5];
-                    aMesh->ChangeElementNodes( f, tmpnodes, nbn );
+                    myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4],
+                                                             nodes[1], nodes[3], nodes[5]));
+                    aMesh->RemoveElement(f);
                   }
                 }
                 else {       /////// quadratic quadrangle
@@ -3454,7 +4098,9 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
                     tmpnodes[5] = nodes[3];
                     tmpnodes[6] = nodes[5];
                     tmpnodes[7] = nodes[7];
-                    aMesh->ChangeElementNodes( f, tmpnodes, nbn );
+                    myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
+                                                             nodes[1], nodes[3], nodes[5], nodes[7]));
+                    aMesh->RemoveElement(f);
                   }
                 }
               }
@@ -3464,7 +4110,11 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
                 if ( !f )
                   myLastCreatedElems.Append(aMesh->AddPolygonalFace(polygon_nodes));
                 else if ( nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
+                  {
+                  // TODO problem ChangeElementNodes : not the same number of nodes, not the same type
+                  MESSAGE("ChangeElementNodes");
                   aMesh->ChangeElementNodes( f, nodes, nbn );
+                  }
               }
             }
             while ( srcElements.Length() < myLastCreatedElems.Length() )
@@ -3762,6 +4412,7 @@ SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet &  theElems,
                                   const int           theFlags,
                                   const double        theTolerance)
 {
+  MESSAGE("ExtrusionSweep " << theMakeGroups << " " << theFlags << " " << theTolerance);
   myLastCreatedElems.Clear();
   myLastCreatedNodes.Clear();
 
@@ -3962,6 +4613,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
                                        const gp_Pnt&        theRefPoint,
                                        const bool           theMakeGroups)
 {
+  MESSAGE("ExtrusionAlongTrack");
   myLastCreatedElems.Clear();
   myLastCreatedNodes.Clear();
 
@@ -4020,7 +4672,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
     while ( aItN->more() ) {
       const SMDS_MeshNode* pNode = aItN->next();
       const SMDS_EdgePosition* pEPos =
-        static_cast<const SMDS_EdgePosition*>( pNode->GetPosition().get() );
+        static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
       double aT = pEPos->GetUParameter();
       aPrms.push_back( aT );
     }
@@ -4064,7 +4716,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
         while ( aItN->more() ) {
           const SMDS_MeshNode* pNode = aItN->next();
           const SMDS_EdgePosition* pEPos =
-            static_cast<const SMDS_EdgePosition*>( pNode->GetPosition().get() );
+            static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
           double aT = pEPos->GetUParameter();
           aPrms.push_back( aT );
         }
@@ -4192,7 +4844,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
       const SMDS_MeshNode* pNode = aItN->next();
       if( pNode==aN1 || pNode==aN2 ) continue;
       const SMDS_EdgePosition* pEPos =
-        static_cast<const SMDS_EdgePosition*>( pNode->GetPosition().get() );
+        static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
       double aT = pEPos->GetUParameter();
       aPrms.push_back( aT );
     }
@@ -4239,7 +4891,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
         while ( aItN->more() ) {
           const SMDS_MeshNode* pNode = aItN->next();
           const SMDS_EdgePosition* pEPos =
-            static_cast<const SMDS_EdgePosition*>( pNode->GetPosition().get() );
+            static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
           double aT = pEPos->GetUParameter();
           aPrms.push_back( aT );
         }
@@ -4365,6 +5017,7 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet&  theElements,
                                    const gp_Pnt& theRefPoint,
                                    const bool theMakeGroups)
 {
+  MESSAGE("MakeExtrElements");
   //cout<<"MakeExtrElements  fullList.size() = "<<fullList.size()<<endl;
   int aNbTP = fullList.size();
   vector<SMESH_MeshEditor_PathPoint> aPPs(aNbTP);
@@ -4520,6 +5173,7 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet&  theElements,
           }
 
           // 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.;
@@ -4641,10 +5295,17 @@ void SMESH_MeshEditor::LinearAngleVariation(const int nbSteps,
 }
 
 
-//=======================================================================
-//function : Transform
-//purpose  :
-//=======================================================================
+//================================================================================
+/*!
+ * \brief Move or copy theElements applying theTrsf to their nodes
+ *  \param theElems - elements to transform, if theElems is empty then apply to all mesh nodes
+ *  \param theTrsf - transformation to apply
+ *  \param theCopy - if true, create translated copies of theElems
+ *  \param theMakeGroups - if true and theCopy, create translated groups
+ *  \param theTargetMesh - mesh to copy translated elements into
+ *  \retval SMESH_MeshEditor::PGroupIDs - list of ids of created groups
+ */
+//================================================================================
 
 SMESH_MeshEditor::PGroupIDs
 SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
@@ -4660,21 +5321,37 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
   string groupPostfix;
   switch ( theTrsf.Form() ) {
   case gp_PntMirror:
+    MESSAGE("gp_PntMirror");
+    needReverse = true;
+    groupPostfix = "mirrored";
+    break;
   case gp_Ax1Mirror:
+    MESSAGE("gp_Ax1Mirror");
+    groupPostfix = "mirrored";
+    break;
   case gp_Ax2Mirror:
+    MESSAGE("gp_Ax2Mirror");
     needReverse = true;
     groupPostfix = "mirrored";
     break;
   case gp_Rotation:
+    MESSAGE("gp_Rotation");
     groupPostfix = "rotated";
     break;
   case gp_Translation:
+    MESSAGE("gp_Translation");
     groupPostfix = "translated";
     break;
   case gp_Scale:
+    MESSAGE("gp_Scale");
+    groupPostfix = "scaled";
+    break;
+  case gp_CompoundTrsf: // different scale by axis
+    MESSAGE("gp_CompoundTrsf");
     groupPostfix = "scaled";
     break;
   default:
+    MESSAGE("default");
     needReverse = false;
     groupPostfix = "transformed";
   }
@@ -4694,9 +5371,29 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
   // source elements for each generated one
   SMESH_SequenceOfElemPtr srcElems, srcNodes;
 
-  // loop on theElems
+  // issue 021015: EDF 1578 SMESH: Free nodes are removed when translating a mesh
+  TIDSortedElemSet orphanNode;
+
+  if ( theElems.empty() ) // transform the whole mesh
+  {
+    // add all elements
+    SMDS_ElemIteratorPtr eIt = aMesh->elementsIterator();
+    while ( eIt->more() ) theElems.insert( eIt->next() );
+    // add orphan nodes
+    SMDS_NodeIteratorPtr nIt = aMesh->nodesIterator();
+    while ( nIt->more() )
+    {
+      const SMDS_MeshNode* node = nIt->next();
+      if ( node->NbInverseElements() == 0)
+        orphanNode.insert( node );
+    }
+  }
+
+  // loop on elements to transform nodes : first orphan nodes then elems
   TIDSortedElemSet::iterator itElem;
-  for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
+  TIDSortedElemSet *elements[] = {&orphanNode, &theElems };
+  for (int i=0; i<2; i++)
+  for ( itElem = elements[i]->begin(); itElem != elements[i]->end(); itElem++ ) {
     const SMDS_MeshElement* elem = *itElem;
     if ( !elem )
       continue;
@@ -4705,8 +5402,8 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
     while ( itN->more() ) {
 
-      // check if a node has been already transformed
       const SMDS_MeshNode* node = cast2Node( itN->next() );
+      // check if a node has been already transformed
       pair<TNodeNodeMap::iterator,bool> n2n_isnew =
         nodeMap.insert( make_pair ( node, node ));
       if ( !n2n_isnew.second )
@@ -4756,7 +5453,7 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
     theElems.insert( *invElemIt );
 
   // replicate or reverse elements
-
+  // TODO revoir ordre reverse vtk
   enum {
     REV_TETRA   = 0,  //  = nbNodes - 4
     REV_PYRAMID = 1,  //  = nbNodes - 4
@@ -4824,8 +5521,8 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
       case SMDSAbs_Volume:
         {
           // ATTENTION: Reversing is not yet done!!!
-          const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
-            dynamic_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
+          const SMDS_VtkVolume* aPolyedre =
+            dynamic_cast<const SMDS_VtkVolume*>( elem );
           if (!aPolyedre) {
             MESSAGE("Warning: bad volumic element");
             continue;
@@ -4872,12 +5569,12 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
 
     // Regular elements
     int* i = index[ FORWARD ];
-    if ( needReverse && nbNodes > 2) // reverse mirrored faces and volumes
+    if ( needReverse && nbNodes > 2) {// reverse mirrored faces and volumes
       if ( elemType == SMDSAbs_Face )
         i = index[ REV_FACE ];
       else
         i = index[ nbNodes - 4 ];
-
+    }
     if(elem->IsQuadratic()) {
       static int anIds[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19};
       i = anIds;
@@ -4936,10 +5633,8 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
       }
     }
     else if ( theCopy ) {
-      if ( SMDS_MeshElement* copy = AddElement( nodes, elem->GetType(), elem->IsPoly() )) {
-        myLastCreatedElems.Append( copy );
+      if ( AddElement( nodes, elem->GetType(), elem->IsPoly() ))
         srcElems.Append( elem );
-      }
     }
     else {
       // reverse element as it was reversed by transformation
@@ -4957,6 +5652,332 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
   return newGroupIDs;
 }
 
+
+////=======================================================================
+////function : Scale
+////purpose  :
+////=======================================================================
+//
+//SMESH_MeshEditor::PGroupIDs
+//SMESH_MeshEditor::Scale (TIDSortedElemSet & theElems,
+//                         const gp_Pnt&            thePoint,
+//                         const std::list<double>& theScaleFact,
+//                         const bool         theCopy,
+//                         const bool         theMakeGroups,
+//                         SMESH_Mesh*        theTargetMesh)
+//{
+//  MESSAGE("Scale");
+//  myLastCreatedElems.Clear();
+//  myLastCreatedNodes.Clear();
+//
+//  SMESH_MeshEditor targetMeshEditor( theTargetMesh );
+//  SMESHDS_Mesh* aTgtMesh = theTargetMesh ? theTargetMesh->GetMeshDS() : 0;
+//  SMESHDS_Mesh* aMesh    = GetMeshDS();
+//
+//  double scaleX=1.0, scaleY=1.0, scaleZ=1.0;
+//  std::list<double>::const_iterator itS = theScaleFact.begin();
+//  scaleX = (*itS);
+//  if(theScaleFact.size()==1) {
+//    scaleY = (*itS);
+//    scaleZ= (*itS);
+//  }
+//  if(theScaleFact.size()==2) {
+//    itS++;
+//    scaleY = (*itS);
+//    scaleZ= (*itS);
+//  }
+//  if(theScaleFact.size()>2) {
+//    itS++;
+//    scaleY = (*itS);
+//    itS++;
+//    scaleZ= (*itS);
+//  }
+//
+//  // map old node to new one
+//  TNodeNodeMap nodeMap;
+//
+//  // elements sharing moved nodes; those of them which have all
+//  // nodes mirrored but are not in theElems are to be reversed
+//  TIDSortedElemSet inverseElemSet;
+//
+//  // source elements for each generated one
+//  SMESH_SequenceOfElemPtr srcElems, srcNodes;
+//
+//  // loop on theElems
+//  TIDSortedElemSet::iterator itElem;
+//  for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
+//    const SMDS_MeshElement* elem = *itElem;
+//    if ( !elem )
+//      continue;
+//
+//    // loop on elem nodes
+//    SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+//    while ( itN->more() ) {
+//
+//      // check if a node has been already transformed
+//      const SMDS_MeshNode* node = cast2Node( itN->next() );
+//      pair<TNodeNodeMap::iterator,bool> n2n_isnew =
+//        nodeMap.insert( make_pair ( node, node ));
+//      if ( !n2n_isnew.second )
+//        continue;
+//
+//      //double coord[3];
+//      //coord[0] = node->X();
+//      //coord[1] = node->Y();
+//      //coord[2] = node->Z();
+//      //theTrsf.Transforms( coord[0], coord[1], coord[2] );
+//      double dx = (node->X() - thePoint.X()) * scaleX;
+//      double dy = (node->Y() - thePoint.Y()) * scaleY;
+//      double dz = (node->Z() - thePoint.Z()) * scaleZ;
+//      if ( theTargetMesh ) {
+//        //const SMDS_MeshNode * newNode = aTgtMesh->AddNode( coord[0], coord[1], coord[2] );
+//        const SMDS_MeshNode * newNode =
+//          aTgtMesh->AddNode( thePoint.X()+dx, thePoint.Y()+dy, thePoint.Z()+dz );
+//        n2n_isnew.first->second = newNode;
+//        myLastCreatedNodes.Append(newNode);
+//        srcNodes.Append( node );
+//      }
+//      else if ( theCopy ) {
+//        //const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
+//        const SMDS_MeshNode * newNode =
+//          aMesh->AddNode( thePoint.X()+dx, thePoint.Y()+dy, thePoint.Z()+dz );
+//        n2n_isnew.first->second = newNode;
+//        myLastCreatedNodes.Append(newNode);
+//        srcNodes.Append( node );
+//      }
+//      else {
+//        //aMesh->MoveNode( node, coord[0], coord[1], coord[2] );
+//        aMesh->MoveNode( node, thePoint.X()+dx, thePoint.Y()+dy, thePoint.Z()+dz );
+//        // node position on shape becomes invalid
+//        const_cast< SMDS_MeshNode* > ( node )->SetPosition
+//          ( SMDS_SpacePosition::originSpacePosition() );
+//      }
+//
+//      // keep inverse elements
+//      //if ( !theCopy && !theTargetMesh && needReverse ) {
+//      //  SMDS_ElemIteratorPtr invElemIt = node->GetInverseElementIterator();
+//      //  while ( invElemIt->more() ) {
+//      //    const SMDS_MeshElement* iel = invElemIt->next();
+//      //    inverseElemSet.insert( iel );
+//      //  }
+//      //}
+//    }
+//  }
+//
+//  // either create new elements or reverse mirrored ones
+//  //if ( !theCopy && !needReverse && !theTargetMesh )
+//  if ( !theCopy && !theTargetMesh )
+//    return PGroupIDs();
+//
+//  TIDSortedElemSet::iterator invElemIt = inverseElemSet.begin();
+//  for ( ; invElemIt != inverseElemSet.end(); invElemIt++ )
+//    theElems.insert( *invElemIt );
+//
+//  // replicate or reverse elements
+//
+//  enum {
+//    REV_TETRA   = 0,  //  = nbNodes - 4
+//    REV_PYRAMID = 1,  //  = nbNodes - 4
+//    REV_PENTA   = 2,  //  = nbNodes - 4
+//    REV_FACE    = 3,
+//    REV_HEXA    = 4,  //  = nbNodes - 4
+//    FORWARD     = 5
+//  };
+//  int index[][8] = {
+//    { 2, 1, 0, 3, 4, 0, 0, 0 },  // REV_TETRA
+//    { 2, 1, 0, 3, 4, 0, 0, 0 },  // REV_PYRAMID
+//    { 2, 1, 0, 5, 4, 3, 0, 0 },  // REV_PENTA
+//    { 2, 1, 0, 3, 0, 0, 0, 0 },  // REV_FACE
+//    { 2, 1, 0, 3, 6, 5, 4, 7 },  // REV_HEXA
+//    { 0, 1, 2, 3, 4, 5, 6, 7 }   // FORWARD
+//  };
+//
+//  for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
+//  {
+//    const SMDS_MeshElement* elem = *itElem;
+//    if ( !elem || elem->GetType() == SMDSAbs_Node )
+//      continue;
+//
+//    int nbNodes = elem->NbNodes();
+//    int elemType = elem->GetType();
+//
+//    if (elem->IsPoly()) {
+//      // Polygon or Polyhedral Volume
+//      switch ( elemType ) {
+//      case SMDSAbs_Face:
+//        {
+//          vector<const SMDS_MeshNode*> poly_nodes (nbNodes);
+//          int iNode = 0;
+//          SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+//          while (itN->more()) {
+//            const SMDS_MeshNode* node =
+//              static_cast<const SMDS_MeshNode*>(itN->next());
+//            TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
+//            if (nodeMapIt == nodeMap.end())
+//              break; // not all nodes transformed
+//            //if (needReverse) {
+//            //  // reverse mirrored faces and volumes
+//            //  poly_nodes[nbNodes - iNode - 1] = (*nodeMapIt).second;
+//            //} else {
+//            poly_nodes[iNode] = (*nodeMapIt).second;
+//            //}
+//            iNode++;
+//          }
+//          if ( iNode != nbNodes )
+//            continue; // not all nodes transformed
+//
+//          if ( theTargetMesh ) {
+//            myLastCreatedElems.Append(aTgtMesh->AddPolygonalFace(poly_nodes));
+//            srcElems.Append( elem );
+//          }
+//          else if ( theCopy ) {
+//            myLastCreatedElems.Append(aMesh->AddPolygonalFace(poly_nodes));
+//            srcElems.Append( elem );
+//          }
+//          else {
+//            aMesh->ChangePolygonNodes(elem, poly_nodes);
+//          }
+//        }
+//        break;
+//      case SMDSAbs_Volume:
+//        {
+//          // ATTENTION: Reversing is not yet done!!!
+//          const SMDS_VtkVolume* aPolyedre =
+//            dynamic_cast<const SMDS_VtkVolume*>( elem );
+//          if (!aPolyedre) {
+//            MESSAGE("Warning: bad volumic element");
+//            continue;
+//          }
+//
+//          vector<const SMDS_MeshNode*> poly_nodes;
+//          vector<int> quantities;
+//
+//          bool allTransformed = true;
+//          int nbFaces = aPolyedre->NbFaces();
+//          for (int iface = 1; iface <= nbFaces && allTransformed; iface++) {
+//            int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
+//            for (int inode = 1; inode <= nbFaceNodes && allTransformed; inode++) {
+//              const SMDS_MeshNode* node = aPolyedre->GetFaceNode(iface, inode);
+//              TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
+//              if (nodeMapIt == nodeMap.end()) {
+//                allTransformed = false; // not all nodes transformed
+//              } else {
+//                poly_nodes.push_back((*nodeMapIt).second);
+//              }
+//            }
+//            quantities.push_back(nbFaceNodes);
+//          }
+//          if ( !allTransformed )
+//            continue; // not all nodes transformed
+//
+//          if ( theTargetMesh ) {
+//            myLastCreatedElems.Append(aTgtMesh->AddPolyhedralVolume(poly_nodes, quantities));
+//            srcElems.Append( elem );
+//          }
+//          else if ( theCopy ) {
+//            myLastCreatedElems.Append(aMesh->AddPolyhedralVolume(poly_nodes, quantities));
+//            srcElems.Append( elem );
+//          }
+//          else {
+//            aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
+//          }
+//        }
+//        break;
+//      default:;
+//      }
+//      continue;
+//    }
+//
+//    // Regular elements
+//    int* i = index[ FORWARD ];
+//    //if ( needReverse && nbNodes > 2) // reverse mirrored faces and volumes
+//    //  if ( elemType == SMDSAbs_Face )
+//    //    i = index[ REV_FACE ];
+//    //  else
+//    //    i = index[ nbNodes - 4 ];
+//
+//    if(elem->IsQuadratic()) {
+//      static int anIds[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19};
+//      i = anIds;
+//      //if(needReverse) {
+//      //  if(nbNodes==3) { // quadratic edge
+//      //    static int anIds[] = {1,0,2};
+//      //    i = anIds;
+//      //  }
+//      //  else if(nbNodes==6) { // quadratic triangle
+//      //    static int anIds[] = {0,2,1,5,4,3};
+//      //    i = anIds;
+//      //  }
+//      //  else if(nbNodes==8) { // quadratic quadrangle
+//      //    static int anIds[] = {0,3,2,1,7,6,5,4};
+//      //    i = anIds;
+//      //  }
+//      //  else if(nbNodes==10) { // quadratic tetrahedron of 10 nodes
+//      //    static int anIds[] = {0,2,1,3,6,5,4,7,9,8};
+//      //    i = anIds;
+//      //  }
+//      //  else if(nbNodes==13) { // quadratic pyramid of 13 nodes
+//      //    static int anIds[] = {0,3,2,1,4,8,7,6,5,9,12,11,10};
+//      //    i = anIds;
+//      //  }
+//      //  else if(nbNodes==15) { // quadratic pentahedron with 15 nodes
+//      //    static int anIds[] = {0,2,1,3,5,4,8,7,6,11,10,9,12,14,13};
+//      //    i = anIds;
+//      //  }
+//      //  else { // nbNodes==20 - quadratic hexahedron with 20 nodes
+//      //    static int anIds[] = {0,3,2,1,4,7,6,5,11,10,9,8,15,14,13,12,16,19,18,17};
+//      //    i = anIds;
+//      //  }
+//      //}
+//    }
+//
+//    // find transformed nodes
+//    vector<const SMDS_MeshNode*> nodes(nbNodes);
+//    int iNode = 0;
+//    SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+//    while ( itN->more() ) {
+//      const SMDS_MeshNode* node =
+//        static_cast<const SMDS_MeshNode*>( itN->next() );
+//      TNodeNodeMap::iterator nodeMapIt = nodeMap.find( node );
+//      if ( nodeMapIt == nodeMap.end() )
+//        break; // not all nodes transformed
+//      nodes[ i [ iNode++ ]] = (*nodeMapIt).second;
+//    }
+//    if ( iNode != nbNodes )
+//      continue; // not all nodes transformed
+//
+//    if ( theTargetMesh ) {
+//      if ( SMDS_MeshElement* copy =
+//           targetMeshEditor.AddElement( nodes, elem->GetType(), elem->IsPoly() )) {
+//        myLastCreatedElems.Append( copy );
+//        srcElems.Append( elem );
+//      }
+//    }
+//    else if ( theCopy ) {
+//      if ( SMDS_MeshElement* copy = AddElement( nodes, elem->GetType(), elem->IsPoly() )) {
+//        myLastCreatedElems.Append( copy );
+//        srcElems.Append( elem );
+//      }
+//    }
+//    else {
+//      // reverse element as it was reversed by transformation
+//      if ( nbNodes > 2 )
+//        aMesh->ChangeElementNodes( elem, &nodes[0], nbNodes );
+//    }
+//  }
+//
+//  PGroupIDs newGroupIDs;
+//
+//  if ( theMakeGroups && theCopy ||
+//       theMakeGroups && theTargetMesh ) {
+//    string groupPostfix = "scaled";
+//    newGroupIDs = generateGroups( srcNodes, srcElems, groupPostfix, theTargetMesh );
+//  }
+//
+//  return newGroupIDs;
+//}
+
+
 //=======================================================================
 /*!
  * \brief Create groups of elements made during transformation
@@ -5091,24 +6112,21 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
  */
 //================================================================================
 
-void SMESH_MeshEditor::FindCoincidentNodes (set<const SMDS_MeshNode*> & theNodes,
-                                            const double                theTolerance,
-                                            TListOfListOfNodes &        theGroupsOfNodes)
+void SMESH_MeshEditor::FindCoincidentNodes (TIDSortedNodeSet &   theNodes,
+                                            const double         theTolerance,
+                                            TListOfListOfNodes & theGroupsOfNodes)
 {
   myLastCreatedElems.Clear();
   myLastCreatedNodes.Clear();
 
-  set<const SMDS_MeshNode*> nodes;
   if ( theNodes.empty() )
   { // get all nodes in the mesh
-    SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator();
+    SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator(/*idInceasingOrder=*/true);
     while ( nIt->more() )
-      nodes.insert( nodes.end(),nIt->next());
+      theNodes.insert( theNodes.end(),nIt->next());
   }
-  else
-    nodes=theNodes;
 
-  SMESH_OctreeNode::FindCoincidentNodes ( nodes, &theGroupsOfNodes, theTolerance);
+  SMESH_OctreeNode::FindCoincidentNodes ( theNodes, &theGroupsOfNodes, theTolerance);
 }
 
 
@@ -5128,9 +6146,9 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
   {
     myMesh = ( SMESHDS_Mesh* ) theMesh;
 
-    set<const SMDS_MeshNode*> nodes;
+    TIDSortedNodeSet nodes;
     if ( theMesh ) {
-      SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator();
+      SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator(/*idInceasingOrder=*/true);
       while ( nIt->more() )
         nodes.insert( nodes.end(), nIt->next() );
     }
@@ -5163,9 +6181,8 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
    */
   const SMDS_MeshNode* FindClosestTo( const gp_Pnt& thePnt )
   {
-    SMDS_MeshNode tgtNode( thePnt.X(), thePnt.Y(), thePnt.Z() );
     map<double, const SMDS_MeshNode*> dist2Nodes;
-    myOctreeNode->NodesAround( &tgtNode, dist2Nodes, myHalfLeafSize );
+    myOctreeNode->NodesAround( thePnt.Coord(), dist2Nodes, myHalfLeafSize );
     if ( !dist2Nodes.empty() )
       return dist2Nodes.begin()->second;
     list<const SMDS_MeshNode*> nodes;
@@ -5181,13 +6198,14 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
       list< SMESH_OctreeNode* >::iterator trIt;
       treeList.push_back( myOctreeNode );
 
-      SMDS_MeshNode pointNode( thePnt.X(), thePnt.Y(), thePnt.Z() );
+      gp_XYZ pointNode( thePnt.X(), thePnt.Y(), thePnt.Z() );
+      bool pointInside = myOctreeNode->isInside( pointNode, myHalfLeafSize );
       for ( trIt = treeList.begin(); trIt != treeList.end(); ++trIt)
       {
         SMESH_OctreeNode* tree = *trIt;
         if ( !tree->isLeaf() ) // put children to the queue
         {
-          if ( !tree->isInside( &pointNode, myHalfLeafSize )) continue;
+          if ( pointInside && !tree->isInside( pointNode, myHalfLeafSize )) continue;
           SMESH_OctreeNodeIteratorPtr cIt = tree->GetChildrenIterator();
           while ( cIt->more() )
             treeList.push_back( cIt->next() );
@@ -5223,7 +6241,7 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
     const SMDS_MeshNode* closestNode = 0;
     list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
     for ( ; nIt != nodes.end(); ++nIt ) {
-      double sqDist = thePnt.SquareDistance( TNodeXYZ( *nIt ) );
+      double sqDist = thePnt.SquareDistance( SMESH_MeshEditor::TNodeXYZ( *nIt ) );
       if ( minSqDist > sqDist ) {
         closestNode = *nIt;
         minSqDist = sqDist;
@@ -5278,8 +6296,9 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
   {
   public:
 
-    ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType);
+    ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, double tolerance = NodeRadius );
     void getElementsNearPoint( const gp_Pnt& point, TIDSortedElemSet& foundElems);
+    void getElementsNearLine ( const gp_Ax1& line, TIDSortedElemSet& foundElems);
     ~ElementBndBoxTree();
 
   protected:
@@ -5293,7 +6312,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
     {
       const SMDS_MeshElement* _element;
       int                     _refCount; // an ElementBox can be included in several tree branches
-      ElementBox(const SMDS_MeshElement* elem);
+      ElementBox(const SMDS_MeshElement* elem, double tolerance);
     };
     vector< ElementBox* > _elements;
   };
@@ -5304,7 +6323,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
    */
   //================================================================================
 
-  ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType)
+  ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, double tolerance)
     :SMESH_Octree( new SMESH_Octree::Limit( MaxLevel, /*minSize=*/0. ))
   {
     int nbElems = mesh.GetMeshInfo().NbElements( elemType );
@@ -5312,7 +6331,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
 
     SMDS_ElemIteratorPtr elemIt = mesh.elementsIterator( elemType );
     while ( elemIt->more() )
-      _elements.push_back( new ElementBox( elemIt->next() ));
+      _elements.push_back( new ElementBox( elemIt->next(),tolerance  ));
 
     if ( _elements.size() > MaxNbElemsInLeaf )
       compute();
@@ -5405,143 +6424,609 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
     }
   }
 
+  //================================================================================
+  /*!
+   * \brief Return elements which can be intersected by the line
+   */
+  //================================================================================
+
+  void ElementBndBoxTree::getElementsNearLine( const gp_Ax1&     line,
+                                               TIDSortedElemSet& foundElems)
+  {
+    if ( level() && getBox().IsOut( line ))
+      return;
+
+    if ( isLeaf() )
+    {
+      for ( int i = 0; i < _elements.size(); ++i )
+        if ( !_elements[i]->IsOut( line ))
+          foundElems.insert( _elements[i]->_element );
+    }
+    else
+    {
+      for (int i = 0; i < 8; i++)
+        ((ElementBndBoxTree*) myChildren[i])->getElementsNearLine( line, foundElems );
+    }
+  }
+
   //================================================================================
   /*!
    * \brief Construct the element box
    */
   //================================================================================
 
-  ElementBndBoxTree::ElementBox::ElementBox(const SMDS_MeshElement* elem)
+  ElementBndBoxTree::ElementBox::ElementBox(const SMDS_MeshElement* elem, double tolerance)
   {
     _element  = elem;
     _refCount = 1;
     SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
     while ( nIt->more() )
-      Add( TNodeXYZ( cast2Node( nIt->next() )));
-    Enlarge( NodeRadius );
+      Add( SMESH_MeshEditor::TNodeXYZ( cast2Node( nIt->next() )));
+    Enlarge( tolerance );
   }
 
 } // namespace
 
 //=======================================================================
 /*!
- * \brief Implementation of search for the elements by point
+ * \brief Implementation of search for the elements by point and
+ *        of classification of point in 2D mesh
  */
 //=======================================================================
 
 struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
 {
-  SMESHDS_Mesh*           _mesh;
-  ElementBndBoxTree*      _ebbTree;
-  SMESH_NodeSearcherImpl* _nodeSearcher;
-  SMDSAbs_ElementType     _elementType;
-
-  SMESH_ElementSearcherImpl( SMESHDS_Mesh& mesh ): _mesh(&mesh),_ebbTree(0),_nodeSearcher(0) {}
+  SMESHDS_Mesh*                _mesh;
+  ElementBndBoxTree*           _ebbTree;
+  SMESH_NodeSearcherImpl*      _nodeSearcher;
+  SMDSAbs_ElementType          _elementType;
+  double                       _tolerance;
+  bool                         _outerFacesFound;
+  set<const SMDS_MeshElement*> _outerFaces; // empty means "no internal faces at all"
+
+  SMESH_ElementSearcherImpl( SMESHDS_Mesh& mesh )
+    : _mesh(&mesh),_ebbTree(0),_nodeSearcher(0), _tolerance(-1), _outerFacesFound(false) {}
   ~SMESH_ElementSearcherImpl()
   {
     if ( _ebbTree )      delete _ebbTree;      _ebbTree      = 0;
     if ( _nodeSearcher ) delete _nodeSearcher; _nodeSearcher = 0;
   }
-
-  /*!
-   * \brief Return elements of given type where the given point is IN or ON.
-   *
-   * 'ALL' type means elements of any type excluding nodes and 0D elements
-   */
-  void FindElementsByPoint(const gp_Pnt&                      point,
-                           SMDSAbs_ElementType                type,
-                           vector< const SMDS_MeshElement* >& foundElements)
+  virtual int FindElementsByPoint(const gp_Pnt&                      point,
+                                  SMDSAbs_ElementType                type,
+                                  vector< const SMDS_MeshElement* >& foundElements);
+  virtual TopAbs_State GetPointState(const gp_Pnt& point);
+
+  void GetElementsNearLine( const gp_Ax1&                      line,
+                            SMDSAbs_ElementType                type,
+                            vector< const SMDS_MeshElement* >& foundElems);
+  double getTolerance();
+  bool getIntersParamOnLine(const gp_Lin& line, const SMDS_MeshElement* face,
+                            const double tolerance, double & param);
+  void findOuterBoundary(const SMDS_MeshElement* anyOuterFace);
+  bool isOuterBoundary(const SMDS_MeshElement* face) const
+  {
+    return _outerFaces.empty() || _outerFaces.count(face);
+  }
+  struct TInters //!< data of intersection of the line and the mesh face (used in GetPointState())
+  {
+    const SMDS_MeshElement* _face;
+    gp_Vec                  _faceNorm;
+    bool                    _coincides; //!< the line lays in face plane
+    TInters(const SMDS_MeshElement* face, const gp_Vec& faceNorm, bool coinc=false)
+      : _face(face), _faceNorm( faceNorm ), _coincides( coinc ) {}
+  };
+  struct TFaceLink //!< link and faces sharing it (used in findOuterBoundary())
   {
-    foundElements.clear();
+    SMESH_TLink      _link;
+    TIDSortedElemSet _faces;
+    TFaceLink( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const SMDS_MeshElement* face)
+      : _link( n1, n2 ), _faces( &face, &face + 1) {}
+  };
+};
+
+ostream& operator<< (ostream& out, const SMESH_ElementSearcherImpl::TInters& i)
+{
+  return out << "TInters(face=" << ( i._face ? i._face->GetID() : 0)
+             << ", _coincides="<<i._coincides << ")";
+}
+
+//=======================================================================
+/*!
+ * \brief define tolerance for search
+ */
+//=======================================================================
 
+double SMESH_ElementSearcherImpl::getTolerance()
+{
+  if ( _tolerance < 0 )
+  {
     const SMDS_MeshInfo& meshInfo = _mesh->GetMeshInfo();
 
-    // -----------------
-    // define tolerance
-    // -----------------
-    double tolerance = 0;
+    _tolerance = 0;
     if ( _nodeSearcher && meshInfo.NbNodes() > 1 )
     {
       double boxSize = _nodeSearcher->getTree()->maxSize();
-      tolerance = 1e-8 * boxSize/* / meshInfo.NbNodes()*/;
+      _tolerance = 1e-8 * boxSize/* / meshInfo.NbNodes()*/;
     }
     else if ( _ebbTree && meshInfo.NbElements() > 0 )
     {
       double boxSize = _ebbTree->maxSize();
-      tolerance = 1e-8 * boxSize/* / meshInfo.NbElements()*/;
+      _tolerance = 1e-8 * boxSize/* / meshInfo.NbElements()*/;
     }
-    if ( tolerance == 0 )
+    if ( _tolerance == 0 )
     {
       // define tolerance by size of a most complex element
       int complexType = SMDSAbs_Volume;
       while ( complexType > SMDSAbs_All &&
               meshInfo.NbElements( SMDSAbs_ElementType( complexType )) < 1 )
         --complexType;
-      if ( complexType == SMDSAbs_All ) return; // empty mesh
-
+      if ( complexType == SMDSAbs_All ) return 0; // empty mesh
       double elemSize;
       if ( complexType == int( SMDSAbs_Node ))
       {
         SMDS_NodeIteratorPtr nodeIt = _mesh->nodesIterator();
         elemSize = 1;
         if ( meshInfo.NbNodes() > 2 )
-          elemSize = TNodeXYZ( nodeIt->next() ).Distance( nodeIt->next() );
+          elemSize = SMESH_MeshEditor::TNodeXYZ( nodeIt->next() ).Distance( nodeIt->next() );
       }
       else
       {
-        const SMDS_MeshElement* elem =
-          _mesh->elementsIterator( SMDSAbs_ElementType( complexType ))->next();
+        SMDS_ElemIteratorPtr elemIt =
+            _mesh->elementsIterator( SMDSAbs_ElementType( complexType ));
+        const SMDS_MeshElement* elem = elemIt->next();
         SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
-        TNodeXYZ n1( cast2Node( nodeIt->next() ));
+        SMESH_MeshEditor::TNodeXYZ n1( cast2Node( nodeIt->next() ));
         while ( nodeIt->more() )
         {
           double dist = n1.Distance( cast2Node( nodeIt->next() ));
           elemSize = max( dist, elemSize );
         }
       }
-      tolerance = 1e-6 * elemSize;
+      _tolerance = 1e-4 * elemSize;
     }
+  }
+  return _tolerance;
+}
+
+//================================================================================
+/*!
+ * \brief Find intersection of the line and an edge of face and return parameter on line
+ */
+//================================================================================
 
-    // =================================================================================
-    if ( type == SMDSAbs_Node || type == SMDSAbs_0DElement )
+bool SMESH_ElementSearcherImpl::getIntersParamOnLine(const gp_Lin&           line,
+                                                     const SMDS_MeshElement* face,
+                                                     const double            tol,
+                                                     double &                param)
+{
+  int nbInts = 0;
+  param = 0;
+
+  GeomAPI_ExtremaCurveCurve anExtCC;
+  Handle(Geom_Curve) lineCurve = new Geom_Line( line );
+  
+  int nbNodes = face->IsQuadratic() ? face->NbNodes()/2 : face->NbNodes();
+  for ( int i = 0; i < nbNodes && nbInts < 2; ++i )
+  {
+    GC_MakeSegment edge( SMESH_MeshEditor::TNodeXYZ( face->GetNode( i )),
+                         SMESH_MeshEditor::TNodeXYZ( face->GetNode( (i+1)%nbNodes) )); 
+    anExtCC.Init( lineCurve, edge);
+    if ( anExtCC.NbExtrema() > 0 && anExtCC.LowerDistance() <= tol)
     {
-      if ( !_nodeSearcher )
-        _nodeSearcher = new SMESH_NodeSearcherImpl( _mesh );
+      Quantity_Parameter pl, pe;
+      anExtCC.LowerDistanceParameters( pl, pe );
+      param += pl;
+      if ( ++nbInts == 2 )
+        break;
+    }
+  }
+  if ( nbInts > 0 ) param /= nbInts;
+  return nbInts > 0;
+}
+//================================================================================
+/*!
+ * \brief Find all faces belonging to the outer boundary of mesh
+ */
+//================================================================================
+
+void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerFace)
+{
+  if ( _outerFacesFound ) return;
+
+  // Collect all outer faces by passing from one outer face to another via their links
+  // and BTW find out if there are internal faces at all.
 
-      const SMDS_MeshNode* closeNode = _nodeSearcher->FindClosestTo( point );
-      if ( !closeNode ) return;
+  // checked links and links where outer boundary meets internal one
+  set< SMESH_TLink > visitedLinks, seamLinks;
 
-      if ( point.Distance( TNodeXYZ( closeNode )) > tolerance )
-        return; // to far from any node
+  // links to treat with already visited faces sharing them
+  list < TFaceLink > startLinks;
 
-      if ( type == SMDSAbs_Node )
+  // load startLinks with the first outerFace
+  startLinks.push_back( TFaceLink( outerFace->GetNode(0), outerFace->GetNode(1), outerFace));
+  _outerFaces.insert( outerFace );
+
+  TIDSortedElemSet emptySet;
+  while ( !startLinks.empty() )
+  {
+    const SMESH_TLink& link  = startLinks.front()._link;
+    TIDSortedElemSet&  faces = startLinks.front()._faces;
+
+    outerFace = *faces.begin();
+    // find other faces sharing the link
+    const SMDS_MeshElement* f;
+    while (( f = SMESH_MeshEditor::FindFaceInSet(link.node1(), link.node2(), emptySet, faces )))
+      faces.insert( f );
+
+    // select another outer face among the found 
+    const SMDS_MeshElement* outerFace2 = 0;
+    if ( faces.size() == 2 )
+    {
+      outerFace2 = (outerFace == *faces.begin() ? *faces.rbegin() : *faces.begin());
+    }
+    else if ( faces.size() > 2 )
+    {
+      seamLinks.insert( link );
+
+      // link direction within the outerFace
+      gp_Vec n1n2( SMESH_MeshEditor::TNodeXYZ( link.node1()),
+                   SMESH_MeshEditor::TNodeXYZ( link.node2()));
+      int i1 = outerFace->GetNodeIndex( link.node1() );
+      int i2 = outerFace->GetNodeIndex( link.node2() );
+      bool rev = ( abs(i2-i1) == 1 ? i1 > i2 : i2 > i1 );
+      if ( rev ) n1n2.Reverse();
+      // outerFace normal
+      gp_XYZ ofNorm, fNorm;
+      if ( SMESH_Algo::FaceNormal( outerFace, ofNorm, /*normalized=*/false ))
+      {
+        // direction from the link inside outerFace
+        gp_Vec dirInOF = gp_Vec( ofNorm ) ^ n1n2;
+        // sort all other faces by angle with the dirInOF
+        map< double, const SMDS_MeshElement* > angle2Face;
+        set< const SMDS_MeshElement*, TIDCompare >::const_iterator face = faces.begin();
+        for ( ; face != faces.end(); ++face )
+        {
+          if ( !SMESH_Algo::FaceNormal( *face, fNorm, /*normalized=*/false ))
+            continue;
+          gp_Vec dirInF = gp_Vec( fNorm ) ^ n1n2;
+          double angle = dirInOF.AngleWithRef( dirInF, n1n2 );
+          if ( angle < 0 ) angle += 2*PI;
+          angle2Face.insert( make_pair( angle, *face ));
+        }
+        if ( !angle2Face.empty() )
+          outerFace2 = angle2Face.begin()->second;
+      }
+    }
+    // store the found outer face and add its links to continue seaching from
+    if ( outerFace2 )
+    {
+      _outerFaces.insert( outerFace );
+      int nbNodes = outerFace2->NbNodes()/( outerFace2->IsQuadratic() ? 2 : 1 );
+      for ( int i = 0; i < nbNodes; ++i )
       {
-        foundElements.push_back( closeNode );
+        SMESH_TLink link2( outerFace2->GetNode(i), outerFace2->GetNode((i+1)%nbNodes));
+        if ( visitedLinks.insert( link2 ).second )
+          startLinks.push_back( TFaceLink( link2.node1(), link2.node2(), outerFace2 ));
       }
-      else
+    }
+    startLinks.pop_front();
+  }
+  _outerFacesFound = true;
+
+  if ( !seamLinks.empty() )
+  {
+    // There are internal boundaries touching the outher one,
+    // find all faces of internal boundaries in order to find
+    // faces of boundaries of holes, if any.
+    
+  }
+  else
+  {
+    _outerFaces.clear();
+  }
+}
+
+//=======================================================================
+/*!
+ * \brief Find elements of given type where the given point is IN or ON.
+ *        Returns nb of found elements and elements them-selves.
+ *
+ * 'ALL' type means elements of any type excluding nodes and 0D elements
+ */
+//=======================================================================
+
+int SMESH_ElementSearcherImpl::
+FindElementsByPoint(const gp_Pnt&                      point,
+                    SMDSAbs_ElementType                type,
+                    vector< const SMDS_MeshElement* >& foundElements)
+{
+  foundElements.clear();
+
+  double tolerance = getTolerance();
+
+  // =================================================================================
+  if ( type == SMDSAbs_Node || type == SMDSAbs_0DElement )
+  {
+    if ( !_nodeSearcher )
+      _nodeSearcher = new SMESH_NodeSearcherImpl( _mesh );
+
+    const SMDS_MeshNode* closeNode = _nodeSearcher->FindClosestTo( point );
+    if ( !closeNode ) return foundElements.size();
+
+    if ( point.Distance( SMESH_MeshEditor::TNodeXYZ( closeNode )) > tolerance )
+      return foundElements.size(); // to far from any node
+
+    if ( type == SMDSAbs_Node )
+    {
+      foundElements.push_back( closeNode );
+    }
+    else
+    {
+      SMDS_ElemIteratorPtr elemIt = closeNode->GetInverseElementIterator( SMDSAbs_0DElement );
+      while ( elemIt->more() )
+        foundElements.push_back( elemIt->next() );
+    }
+  }
+  // =================================================================================
+  else // elements more complex than 0D
+  {
+    if ( !_ebbTree || _elementType != type )
+    {
+      if ( _ebbTree ) delete _ebbTree;
+      _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, tolerance );
+    }
+    TIDSortedElemSet suspectElems;
+    _ebbTree->getElementsNearPoint( point, suspectElems );
+    TIDSortedElemSet::iterator elem = suspectElems.begin();
+    for ( ; elem != suspectElems.end(); ++elem )
+      if ( !SMESH_MeshEditor::isOut( *elem, point, tolerance ))
+        foundElements.push_back( *elem );
+  }
+  return foundElements.size();
+}
+
+//================================================================================
+/*!
+ * \brief Classify the given point in the closed 2D mesh
+ */
+//================================================================================
+
+TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
+{
+  double tolerance = getTolerance();
+  if ( !_ebbTree || _elementType != SMDSAbs_Face )
+  {
+    if ( _ebbTree ) delete _ebbTree;
+    _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = SMDSAbs_Face );
+  }
+  // Algo: analyse transition of a line starting at the point through mesh boundary;
+  // try three lines parallel to axis of the coordinate system and perform rough
+  // analysis. If solution is not clear perform thorough analysis.
+
+  const int nbAxes = 3;
+  gp_Dir axisDir[ nbAxes ] = { gp::DX(), gp::DY(), gp::DZ() };
+  map< double, TInters >   paramOnLine2TInters[ nbAxes ];
+  list< TInters > tangentInters[ nbAxes ]; // of faces whose plane includes the line
+  multimap< int, int > nbInt2Axis; // to find the simplest case
+  for ( int axis = 0; axis < nbAxes; ++axis )
+  {
+    gp_Ax1 lineAxis( point, axisDir[axis]);
+    gp_Lin line    ( lineAxis );
+
+    TIDSortedElemSet suspectFaces; // faces possibly intersecting the line
+    _ebbTree->getElementsNearLine( lineAxis, suspectFaces );
+
+    // Intersect faces with the line
+
+    map< double, TInters > & u2inters = paramOnLine2TInters[ axis ];
+    TIDSortedElemSet::iterator face = suspectFaces.begin();
+    for ( ; face != suspectFaces.end(); ++face )
+    {
+      // get face plane
+      gp_XYZ fNorm;
+      if ( !SMESH_Algo::FaceNormal( *face, fNorm, /*normalized=*/false)) continue;
+      gp_Pln facePlane( SMESH_MeshEditor::TNodeXYZ( (*face)->GetNode(0)), fNorm );
+
+      // perform intersection
+      IntAna_IntConicQuad intersection( line, IntAna_Quadric( facePlane ));
+      if ( !intersection.IsDone() )
+        continue;
+      if ( intersection.IsInQuadric() )
+      {
+        tangentInters[ axis ].push_back( TInters( *face, fNorm, true ));
+      }
+      else if ( ! intersection.IsParallel() && intersection.NbPoints() > 0 )
       {
-        SMDS_ElemIteratorPtr elemIt = closeNode->GetInverseElementIterator( SMDSAbs_0DElement );
-        while ( elemIt->more() )
-          foundElements.push_back( elemIt->next() );
+        gp_Pnt intersectionPoint = intersection.Point(1);
+        if ( !SMESH_MeshEditor::isOut( *face, intersectionPoint, tolerance ))
+          u2inters.insert(make_pair( intersection.ParamOnConic(1), TInters( *face, fNorm )));
       }
     }
-    // =================================================================================
-    else // elements more complex than 0D
+    // Analyse intersections roughly
+
+    int nbInter = u2inters.size();
+    if ( nbInter == 0 )
+      return TopAbs_OUT; 
+
+    double f = u2inters.begin()->first, l = u2inters.rbegin()->first;
+    if ( nbInter == 1 ) // not closed mesh
+      return fabs( f ) < tolerance ? TopAbs_ON : TopAbs_UNKNOWN;
+
+    if ( fabs( f ) < tolerance || fabs( l ) < tolerance )
+      return TopAbs_ON;
+
+    if ( (f<0) == (l<0) )
+      return TopAbs_OUT;
+
+    int nbIntBeforePoint = std::distance( u2inters.begin(), u2inters.lower_bound(0));
+    int nbIntAfterPoint  = nbInter - nbIntBeforePoint;
+    if ( nbIntBeforePoint == 1 || nbIntAfterPoint == 1 )
+      return TopAbs_IN;
+
+    nbInt2Axis.insert( make_pair( min( nbIntBeforePoint, nbIntAfterPoint ), axis ));
+
+    if ( _outerFacesFound ) break; // pass to thorough analysis
+
+  } // three attempts - loop on CS axes
+
+  // Analyse intersections thoroughly.
+  // We make two loops maximum, on the first one we only exclude touching intersections,
+  // on the second, if situation is still unclear, we gather and use information on
+  // position of faces (internal or outer). If faces position is already gathered,
+  // we make the second loop right away.
+
+  for ( int hasPositionInfo = _outerFacesFound; hasPositionInfo < 2; ++hasPositionInfo )
+  {
+    multimap< int, int >::const_iterator nb_axis = nbInt2Axis.begin();
+    for ( ; nb_axis != nbInt2Axis.end(); ++nb_axis )
     {
-      if ( !_ebbTree || _elementType != type )
+      int axis = nb_axis->second;
+      map< double, TInters > & u2inters = paramOnLine2TInters[ axis ];
+
+      gp_Ax1 lineAxis( point, axisDir[axis]);
+      gp_Lin line    ( lineAxis );
+
+      // add tangent intersections to u2inters
+      double param;
+      list< TInters >::const_iterator tgtInt = tangentInters[ axis ].begin();
+      for ( ; tgtInt != tangentInters[ axis ].end(); ++tgtInt )
+        if ( getIntersParamOnLine( line, tgtInt->_face, tolerance, param ))
+          u2inters.insert(make_pair( param, *tgtInt ));
+      tangentInters[ axis ].clear();
+
+      // Count intersections before and after the point excluding touching ones.
+      // If hasPositionInfo we count intersections of outer boundary only
+
+      int nbIntBeforePoint = 0, nbIntAfterPoint = 0;
+      double f = numeric_limits<double>::max(), l = -numeric_limits<double>::max();
+      map< double, TInters >::iterator u_int1 = u2inters.begin(), u_int2 = u_int1;
+      bool ok = ! u_int1->second._coincides;
+      while ( ok && u_int1 != u2inters.end() )
+      {
+        double u = u_int1->first;
+        bool touchingInt = false;
+        if ( ++u_int2 != u2inters.end() )
+        {
+          // skip intersections at the same point (if the line passes through edge or node)
+          int nbSamePnt = 0;
+          while ( u_int2 != u2inters.end() && fabs( u_int2->first - u ) < tolerance )
+          {
+            ++nbSamePnt;
+            ++u_int2;
+          }
+
+          // skip tangent intersections
+          int nbTgt = 0;
+          const SMDS_MeshElement* prevFace = u_int1->second._face;
+          while ( ok && u_int2->second._coincides )
+          {
+            if ( SMESH_Algo::GetCommonNodes(prevFace , u_int2->second._face).empty() )
+              ok = false;
+            else
+            {
+              nbTgt++;
+              u_int2++;
+              ok = ( u_int2 != u2inters.end() );
+            }
+          }
+          if ( !ok ) break;
+
+          // skip intersections at the same point after tangent intersections
+          if ( nbTgt > 0 )
+          {
+            double u2 = u_int2->first;
+            ++u_int2;
+            while ( u_int2 != u2inters.end() && fabs( u_int2->first - u2 ) < tolerance )
+            {
+              ++nbSamePnt;
+              ++u_int2;
+            }
+          }
+          // decide if we skipped a touching intersection
+          if ( nbSamePnt + nbTgt > 0 )
+          {
+            double minDot = numeric_limits<double>::max(), maxDot = -numeric_limits<double>::max();
+            map< double, TInters >::iterator u_int = u_int1;
+            for ( ; u_int != u_int2; ++u_int )
+            {
+              if ( u_int->second._coincides ) continue;
+              double dot = u_int->second._faceNorm * line.Direction();
+              if ( dot > maxDot ) maxDot = dot;
+              if ( dot < minDot ) minDot = dot;
+            }
+            touchingInt = ( minDot*maxDot < 0 );
+          }
+        }
+        if ( !touchingInt )
+        {
+          if ( !hasPositionInfo || isOuterBoundary( u_int1->second._face ))
+          {
+            if ( u < 0 )
+              ++nbIntBeforePoint;
+            else
+              ++nbIntAfterPoint;
+          }
+          if ( u < f ) f = u;
+          if ( u > l ) l = u;
+        }
+
+        u_int1 = u_int2; // to next intersection
+
+      } // loop on intersections with one line
+
+      if ( ok )
       {
-        if ( _ebbTree ) delete _ebbTree;
-        _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type );
+        if ( fabs( f ) < tolerance || fabs( l ) < tolerance )
+          return TopAbs_ON;
+
+        if ( nbIntBeforePoint == 0  || nbIntAfterPoint == 0)
+          return TopAbs_OUT; 
+
+        if ( nbIntBeforePoint + nbIntAfterPoint == 1 ) // not closed mesh
+          return fabs( f ) < tolerance ? TopAbs_ON : TopAbs_UNKNOWN;
+
+        if ( nbIntBeforePoint == 1 || nbIntAfterPoint == 1 )
+          return TopAbs_IN;
+
+        if ( (f<0) == (l<0) )
+          return TopAbs_OUT;
+
+        if ( hasPositionInfo )
+          return nbIntBeforePoint % 2 ? TopAbs_IN : TopAbs_OUT;
       }
-      TIDSortedElemSet suspectElems;
-      _ebbTree->getElementsNearPoint( point, suspectElems );
-      TIDSortedElemSet::iterator elem = suspectElems.begin();
-      for ( ; elem != suspectElems.end(); ++elem )
-        if ( !SMESH_MeshEditor::isOut( *elem, point, tolerance ))
-          foundElements.push_back( *elem );
+    } // loop on intersections of the tree lines - thorough analysis
+
+    if ( !hasPositionInfo )
+    {
+      // gather info on faces position - is face in the outer boundary or not
+      map< double, TInters > & u2inters = paramOnLine2TInters[ 0 ];
+      findOuterBoundary( u2inters.begin()->second._face );
     }
+
+  } // two attempts - with and w/o faces position info in the mesh
+
+  return TopAbs_UNKNOWN;
+}
+
+//=======================================================================
+/*!
+ * \brief Return elements possibly intersecting the line
+ */
+//=======================================================================
+
+void SMESH_ElementSearcherImpl::GetElementsNearLine( const gp_Ax1&                      line,
+                                                     SMDSAbs_ElementType                type,
+                                                     vector< const SMDS_MeshElement* >& foundElems)
+{
+  if ( !_ebbTree || _elementType != type )
+  {
+    if ( _ebbTree ) delete _ebbTree;
+    _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type );
   }
-}; // struct SMESH_ElementSearcherImpl
+  TIDSortedElemSet suspectFaces; // elements possibly intersecting the line
+  _ebbTree->getElementsNearLine( line, suspectFaces );
+  foundElems.assign( suspectFaces.begin(), suspectFaces.end());
+}
 
 //=======================================================================
 /*!
@@ -5570,86 +7055,140 @@ bool SMESH_MeshEditor::isOut( const SMDS_MeshElement* element, const gp_Pnt& poi
   // get ordered nodes
 
   vector< gp_XYZ > xyz;
+  vector<const SMDS_MeshNode*> nodeList;
 
   SMDS_ElemIteratorPtr nodeIt = element->nodesIterator();
-  if ( element->IsQuadratic() )
-    if (const SMDS_QuadraticFaceOfNodes* f=dynamic_cast<const SMDS_QuadraticFaceOfNodes*>(element))
+  if ( element->IsQuadratic() ) {
+    if (const SMDS_VtkFace* f=dynamic_cast<const SMDS_VtkFace*>(element))
       nodeIt = f->interlacedNodesElemIterator();
-    else if (const SMDS_QuadraticEdge*  e =dynamic_cast<const SMDS_QuadraticEdge*>(element))
+    else if (const SMDS_VtkEdge*  e =dynamic_cast<const SMDS_VtkEdge*>(element))
       nodeIt = e->interlacedNodesElemIterator();
-
+  }
   while ( nodeIt->more() )
-    xyz.push_back( TNodeXYZ( cast2Node( nodeIt->next() )));
+    {
+      const SMDS_MeshNode* node = cast2Node( nodeIt->next() );
+      xyz.push_back( TNodeXYZ(node) );
+      nodeList.push_back(node);
+    }
+
+  int i, nbNodes = element->NbNodes();
 
   if ( element->GetType() == SMDSAbs_Face ) // --------------------------------------------------
   {
-    // gravity center
-    gp_XYZ gc(0,0,0);
-    gc = accumulate( xyz.begin(), xyz.end(), gc );
-    gc /= element->NbNodes();
-
-    // compute face normal using gc
-    gp_Vec normal(0,0,0);
+    // compute face normal
+    gp_Vec faceNorm(0,0,0);
     xyz.push_back( xyz.front() );
-    for ( int i = 0; i < element->NbNodes(); ++i )
+    nodeList.push_back( nodeList.front() );
+    for ( i = 0; i < nbNodes; ++i )
+    {
+      gp_Vec edge1( xyz[i+1], xyz[i]);
+      gp_Vec edge2( xyz[i+1], xyz[(i+2)%nbNodes] );
+      faceNorm += edge1 ^ edge2;
+    }
+    double normSize = faceNorm.Magnitude();
+    if ( normSize <= tol )
     {
-      gp_Vec edge( xyz[i], xyz[i+1]);
-      gp_Vec n2gc( xyz[i], gc );
-      normal += edge ^ n2gc;
+      // degenerated face: point is out if it is out of all face edges
+      for ( i = 0; i < nbNodes; ++i )
+      {
+        SMDS_LinearEdge edge( nodeList[i], nodeList[i+1] );
+        if ( !isOut( &edge, point, tol ))
+          return false;
+      }
+      return true;
     }
-    double faceDoubleArea = normal.Magnitude();
-    if ( faceDoubleArea <= numeric_limits<double>::min() )
-      return true; // invalid face
-    normal /= faceDoubleArea;
+    faceNorm /= normSize;
 
     // check if the point lays on face plane
     gp_Vec n2p( xyz[0], point );
-    if ( fabs( n2p * normal ) > tol )
+    if ( fabs( n2p * faceNorm ) > tol )
       return true; // not on face plane
 
-    // check if point is out of face boundary
-    int i, out = false;
-    for ( i = 0; !out && i < element->NbNodes(); ++i )
+    // check if point is out of face boundary:
+    // define it by closest transition of a ray point->infinity through face boundary
+    // on the face plane.
+    // First, find normal of a plane perpendicular to face plane, to be used as a cutting tool
+    // to find intersections of the ray with the boundary.
+    gp_Vec ray = n2p;
+    gp_Vec plnNorm = ray ^ faceNorm;
+    normSize = plnNorm.Magnitude();
+    if ( normSize <= tol ) return false; // point coincides with the first node
+    plnNorm /= normSize;
+    // for each node of the face, compute its signed distance to the plane
+    vector<double> dist( nbNodes + 1);
+    for ( i = 0; i < nbNodes; ++i )
     {
-      gp_Vec edge( xyz[i], xyz[i+1]);
-      gp_Vec n2p ( xyz[i], point );
-      gp_Vec cross = edge ^ n2p;
-      out = ( cross * normal < -tol );
-    }
-    if ( out && element->IsPoly() )
+      gp_Vec n2p( xyz[i], point );
+      dist[i] = n2p * plnNorm;
+    }
+    dist.back() = dist.front();
+    // find the closest intersection
+    int    iClosest = -1;
+    double rClosest, distClosest = 1e100;;
+    gp_Pnt pClosest;
+    for ( i = 0; i < nbNodes; ++i )
     {
-      // define point position by the closest edge
-      double minDist = numeric_limits<double>::max();
-      int iMinDist;
-      for ( i = 0; i < element->NbNodes(); ++i )
+      double r;
+      if ( fabs( dist[i]) < tol )
+        r = 0.;
+      else if ( fabs( dist[i+1]) < tol )
+        r = 1.;
+      else if ( dist[i] * dist[i+1] < 0 )
+        r = dist[i] / ( dist[i] - dist[i+1] );
+      else
+        continue; // no intersection
+      gp_Pnt pInt = xyz[i] * (1.-r) + xyz[i+1] * r;
+      gp_Vec p2int ( point, pInt);
+      if ( p2int * ray > -tol ) // right half-space
       {
-        gp_Vec edge( xyz[i], xyz[i+1]);
-        gp_Vec n1p ( xyz[i], point);
-        double dist = ( edge ^ n1p ).Magnitude() / edge.Magnitude();
-        if ( dist < minDist )
-          iMinDist = i;
+        double intDist = p2int.SquareMagnitude();
+        if ( intDist < distClosest )
+        {
+          iClosest = i;
+          rClosest = r;
+          pClosest = pInt;
+          distClosest = intDist;
+        }
       }
-      gp_Vec edge( xyz[iMinDist], xyz[iMinDist+1]);
-      gp_Vec n2p ( xyz[iMinDist], point );
-      gp_Vec cross = edge ^ n2p;
-      out = ( cross * normal < -tol );
     }
-    return out;
+    if ( iClosest < 0 )
+      return true; // no intesections - out
+
+    // analyse transition
+    gp_Vec edge( xyz[iClosest], xyz[iClosest+1] );
+    gp_Vec edgeNorm = -( edge ^ faceNorm ); // normal to intersected edge pointing out of face
+    gp_Vec p2int ( point, pClosest );
+    bool out = (edgeNorm * p2int) < -tol;
+    if ( rClosest > 0. && rClosest < 1. ) // not node intersection
+      return out;
+
+    // ray pass through a face node; analyze transition through an adjacent edge
+    gp_Pnt p1 = xyz[ (rClosest == 0.) ? ((iClosest+nbNodes-1) % nbNodes) : (iClosest+1) ];
+    gp_Pnt p2 = xyz[ (rClosest == 0.) ? iClosest : ((iClosest+2) % nbNodes) ];
+    gp_Vec edgeAdjacent( p1, p2 );
+    gp_Vec edgeNorm2 = -( edgeAdjacent ^ faceNorm );
+    bool out2 = (edgeNorm2 * p2int) < -tol;
+
+    bool covexCorner = ( edgeNorm * edgeAdjacent * (rClosest==1. ? 1. : -1.)) < 0;
+    return covexCorner ? (out || out2) : (out && out2);
   }
   if ( element->GetType() == SMDSAbs_Edge ) // --------------------------------------------------
   {
-    for ( int i = 1; i < element->NbNodes(); ++i )
+    // point is out of edge if it is NOT ON any straight part of edge
+    // (we consider quadratic edge as being composed of two straight parts)
+    for ( i = 1; i < nbNodes; ++i )
     {
       gp_Vec edge( xyz[i-1], xyz[i]);
       gp_Vec n1p ( xyz[i-1], point);
       double dist = ( edge ^ n1p ).Magnitude() / edge.Magnitude();
       if ( dist > tol )
-        return true;
+        continue;
       gp_Vec n2p( xyz[i], point );
       if ( fabs( edge.Magnitude() - n1p.Magnitude() - n2p.Magnitude()) > tol )
-        return true;
+        continue;
+      return false; // point is ON this part
     }
-    return false;
+    return true;
   }
   // Node or 0D element -------------------------------------------------------------------------
   {
@@ -5750,6 +7289,7 @@ int SMESH_MeshEditor::SimplifyFace (const vector<const SMDS_MeshNode *> faceNode
 
 void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
 {
+  MESSAGE("MergeNodes");
   myLastCreatedElems.Clear();
   myLastCreatedNodes.Clear();
 
@@ -5766,10 +7306,12 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
     list<const SMDS_MeshNode*>& nodes = *grIt;
     list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
     const SMDS_MeshNode* nToKeep = *nIt;
+    //MESSAGE("node to keep " << nToKeep->GetID());
     for ( ++nIt; nIt != nodes.end(); nIt++ ) {
       const SMDS_MeshNode* nToRemove = *nIt;
       nodeNodeMap.insert( TNodeNodeMap::value_type( nToRemove, nToKeep ));
       if ( nToRemove != nToKeep ) {
+        //MESSAGE("  node to remove " << nToRemove->GetID());
         rmNodeIds.push_back( nToRemove->GetID() );
         AddToSameGroups( nToKeep, nToRemove, aMesh );
       }
@@ -5786,6 +7328,7 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
   set<const SMDS_MeshElement*>::iterator eIt = elems.begin();
   for ( ; eIt != elems.end(); eIt++ ) {
     const SMDS_MeshElement* elem = *eIt;
+    //MESSAGE(" ---- inverse elem on node to remove " << elem->GetID());
     int nbNodes = elem->NbNodes();
     int aShapeId = FindShape( elem );
 
@@ -5835,6 +7378,7 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
 
     bool isOk = true;
     int nbUniqueNodes = nodeSet.size();
+    //MESSAGE("nbNodes nbUniqueNodes " << nbNodes << " " << nbUniqueNodes);
     if ( nbNodes != nbUniqueNodes ) { // some nodes stick
       // Polygons and Polyhedral volumes
       if (elem->IsPoly()) {
@@ -5850,10 +7394,9 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
           vector<const SMDS_MeshNode *> polygons_nodes;
           vector<int> quantities;
           int nbNew = SimplifyFace(face_nodes, polygons_nodes, quantities);
-
           if (nbNew > 0) {
             inode = 0;
-            for (int iface = 0; iface < nbNew - 1; iface++) {
+            for (int iface = 0; iface < nbNew; iface++) {
               int nbNodes = quantities[iface];
               vector<const SMDS_MeshNode *> poly_nodes (nbNodes);
               for (int ii = 0; ii < nbNodes; ii++, inode++) {
@@ -5864,7 +7407,19 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
               if (aShapeId)
                 aMesh->SetMeshElementOnShape(newElem, aShapeId);
             }
-            aMesh->ChangeElementNodes(elem, &polygons_nodes[inode], quantities[nbNew - 1]);
+
+            MESSAGE("ChangeElementNodes MergeNodes Polygon");
+            //aMesh->ChangeElementNodes(elem, &polygons_nodes[inode], quantities[nbNew - 1]);
+            vector<const SMDS_MeshNode *> polynodes(polygons_nodes.begin()+inode,polygons_nodes.end());
+            int quid =0;
+            if (nbNew > 0) quid = nbNew - 1;
+            vector<int> newquant(quantities.begin()+quid, quantities.end());
+            const SMDS_MeshElement* newElem = 0;
+            newElem = aMesh->AddPolyhedralVolume(polynodes, newquant);
+            myLastCreatedElems.Append(newElem);
+            if ( aShapeId && newElem )
+              aMesh->SetMeshElementOnShape( newElem, aShapeId );
+            rmElemIds.push_back(elem->GetID());
           }
           else {
             rmElemIds.push_back(elem->GetID());
@@ -5877,9 +7432,9 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
             rmElemIds.push_back(elem->GetID());
           }
           else {
-            // each face has to be analized in order to check volume validity
-            const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
-              static_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
+            // each face has to be analyzed in order to check volume validity
+            const SMDS_VtkVolume* aPolyedre =
+              dynamic_cast<const SMDS_VtkVolume*>( elem );
             if (aPolyedre) {
               int nbFaces = aPolyedre->NbFaces();
 
@@ -5907,10 +7462,16 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
               }
 
               if (quantities.size() > 3)
-                aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
-              else
-                rmElemIds.push_back(elem->GetID());
-
+                {
+                  MESSAGE("ChangeElementNodes MergeNodes Polyhedron");
+                  //aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
+                  const SMDS_MeshElement* newElem = 0;
+                  newElem = aMesh->AddPolyhedralVolume(poly_nodes, quantities);
+                  myLastCreatedElems.Append(newElem);
+                  if ( aShapeId && newElem )
+                    aMesh->SetMeshElementOnShape( newElem, aShapeId );
+                  rmElemIds.push_back(elem->GetID());
+                }
             }
             else {
               rmElemIds.push_back(elem->GetID());
@@ -5924,6 +7485,7 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
       }
 
       // Regular elements
+      // TODO not all the possible cases are solved. Find something more generic?
       switch ( nbNodes ) {
       case 2: ///////////////////////////////////// EDGE
         isOk = false; break;
@@ -5937,6 +7499,7 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
             isOk = false;
           else if ( nbRepl == 2 && iRepl[ 1 ] - iRepl[ 0 ] == 2 )
             isOk = false; // opposite nodes stick
+          //MESSAGE("isOk " << isOk);
         }
         break;
       case 6: ///////////////////////////////////// PENTAHEDRON
@@ -6001,7 +7564,11 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
           //    +---+---+
           //   0    7    3
           isOk = false;
+          if(nbRepl==2) {
+            MESSAGE("nbRepl=2: " << iRepl[0] << " " << iRepl[1]);
+          }
           if(nbRepl==3) {
+            MESSAGE("nbRepl=3: " << iRepl[0] << " " << iRepl[1]  << " " << iRepl[2]);
             nbUniqueNodes = 6;
             if( iRepl[0]==0 && iRepl[1]==1 && iRepl[2]==4 ) {
               uniqueNodes[0] = curNodes[0];
@@ -6076,6 +7643,12 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
               isOk = true;
             }
           }
+          if(nbRepl==4) {
+            MESSAGE("nbRepl=4: " << iRepl[0] << " " << iRepl[1]  << " " << iRepl[2] << " " << iRepl[3]);
+          }
+          if(nbRepl==5) {
+            MESSAGE("nbRepl=5: " << iRepl[0] << " " << iRepl[1]  << " " << iRepl[2] << " " << iRepl[3] << " " << iRepl[4]);
+          }
           break;
         }
         //////////////////////////////////// HEXAHEDRON
@@ -6260,8 +7833,8 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
     if ( isOk ) {
       if (elem->IsPoly() && elem->GetType() == SMDSAbs_Volume) {
         // Change nodes of polyedre
-        const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
-          static_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
+        const SMDS_VtkVolume* aPolyedre =
+          dynamic_cast<const SMDS_VtkVolume*>( elem );
         if (aPolyedre) {
           int nbFaces = aPolyedre->NbFaces();
 
@@ -6286,21 +7859,36 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
         }
       }
       else {
-        // Change regular element or polygon
-        aMesh->ChangeElementNodes( elem, & uniqueNodes[0], nbUniqueNodes );
+        //int elemId = elem->GetID();
+        //MESSAGE("Change regular element or polygon " << elemId);
+        SMDSAbs_ElementType etyp = elem->GetType();
+        uniqueNodes.resize(nbUniqueNodes);
+        SMDS_MeshElement* newElem = 0;
+        if (elem->GetEntityType() == SMDSEntity_Polygon)
+          newElem = this->AddElement(uniqueNodes, etyp, true);
+        else
+          newElem = this->AddElement(uniqueNodes, etyp, false);
+        if (newElem)
+          {
+            myLastCreatedElems.Append(newElem);
+            if ( aShapeId )
+              aMesh->SetMeshElementOnShape( newElem, aShapeId );
+          }
+        aMesh->RemoveElement(elem);
       }
     }
     else {
       // Remove invalid regular element or invalid polygon
+      //MESSAGE("Remove invalid " << elem->GetID());
       rmElemIds.push_back( elem->GetID() );
     }
 
   } // loop on elements
 
-  // Remove equal nodes and bad elements
+  // Remove bad elements, then equal nodes (order important)
 
-  Remove( rmNodeIds, true );
   Remove( rmElemIds, false );
+  Remove( rmNodeIds, true );
 
 }
 
@@ -6447,83 +8035,68 @@ void SMESH_MeshEditor::MergeEqualElements()
 //purpose  : Return a face having linked nodes n1 and n2 and which is
 //           - not in avoidSet,
 //           - in elemSet provided that !elemSet.empty()
+//           i1 and i2 optionally returns indices of n1 and n2
 //=======================================================================
 
 const SMDS_MeshElement*
 SMESH_MeshEditor::FindFaceInSet(const SMDS_MeshNode*    n1,
                                 const SMDS_MeshNode*    n2,
                                 const TIDSortedElemSet& elemSet,
-                                const TIDSortedElemSet& avoidSet)
+                                const TIDSortedElemSet& avoidSet,
+                                int*                    n1ind,
+                                int*                    n2ind)
 
 {
+  int i1, i2;
+  const SMDS_MeshElement* face = 0;
+
   SMDS_ElemIteratorPtr invElemIt = n1->GetInverseElementIterator(SMDSAbs_Face);
-  while ( invElemIt->more() ) { // loop on inverse elements of n1
+  //MESSAGE("n1->GetInverseElementIterator(SMDSAbs_Face) " << invElemIt);
+  while ( invElemIt->more() && !face ) // loop on inverse faces of n1
+  {
+    //MESSAGE("in while ( invElemIt->more() && !face )");
     const SMDS_MeshElement* elem = invElemIt->next();
-    if (avoidSet.find( elem ) != avoidSet.end() )
+    if (avoidSet.count( elem ))
       continue;
-    if ( !elemSet.empty() && elemSet.find( elem ) == elemSet.end())
+    if ( !elemSet.empty() && !elemSet.count( elem ))
       continue;
-    // get face nodes and find index of n1
-    int i1, nbN = elem->NbNodes(), iNode = 0;
-    //const SMDS_MeshNode* faceNodes[ nbN ], *n;
-    vector<const SMDS_MeshNode*> faceNodes( nbN );
-    const SMDS_MeshNode* n;
-    SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
-    while ( nIt->more() ) {
-      faceNodes[ iNode ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
-      if ( faceNodes[ iNode++ ] == n1 )
-        i1 = iNode - 1;
-    }
+    // index of n1
+    i1 = elem->GetNodeIndex( n1 );
     // find a n2 linked to n1
-    if(!elem->IsQuadratic()) {
-      for ( iNode = 0; iNode < 2; iNode++ ) {
-        if ( iNode ) // node before n1
-          n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ];
-        else         // node after n1
-          n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ];
-        if ( n == n2 )
-          return elem;
-      }
-    }
-    else { // analysis for quadratic elements
-      bool IsFind = false;
-      // check using only corner nodes
-      for ( iNode = 0; iNode < 2; iNode++ ) {
-        if ( iNode ) // node before n1
-          n = faceNodes[ i1 == 0 ? nbN/2 - 1 : i1 - 1 ];
-        else         // node after n1
-          n = faceNodes[ i1 + 1 == nbN/2 ? 0 : i1 + 1 ];
-        if ( n == n2 )
-          IsFind = true;
-      }
-      if(IsFind) {
-        return elem;
-      }
-      else {
-        // check using all nodes
-        const SMDS_QuadraticFaceOfNodes* F =
-          static_cast<const SMDS_QuadraticFaceOfNodes*>(elem);
-        // use special nodes iterator
-        iNode = 0;
-        SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
-        while ( anIter->more() ) {
-          faceNodes[iNode] = static_cast<const SMDS_MeshNode*>(anIter->next());
-          if ( faceNodes[ iNode++ ] == n1 )
-            i1 = iNode - 1;
-        }
-        for ( iNode = 0; iNode < 2; iNode++ ) {
-          if ( iNode ) // node before n1
-            n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ];
-          else         // node after n1
-            n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ];
-          if ( n == n2 ) {
-            return elem;
-          }
+    int nbN = elem->IsQuadratic() ? elem->NbNodes()/2 : elem->NbNodes();
+    for ( int di = -1; di < 2 && !face; di += 2 )
+    {
+      i2 = (i1+di+nbN) % nbN;
+      if ( elem->GetNode( i2 ) == n2 )
+        face = elem;
+    }
+    if ( !face && elem->IsQuadratic())
+    {
+      // analysis for quadratic elements using all nodes
+      const SMDS_VtkFace* F =
+        dynamic_cast<const SMDS_VtkFace*>(elem);
+      if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
+      // use special nodes iterator
+      SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
+      const SMDS_MeshNode* prevN = cast2Node( anIter->next() );
+      for ( i1 = -1, i2 = 0; anIter->more() && !face; i1++, i2++ )
+      {
+        const SMDS_MeshNode* n = cast2Node( anIter->next() );
+        if ( n1 == prevN && n2 == n )
+        {
+          face = elem;
+        }
+        else if ( n2 == prevN && n1 == n )
+        {
+          face = elem; swap( i1, i2 );
         }
+        prevN = n;
       }
-    } // end analysis for quadratic elements
+    }
   }
-  return 0;
+  if ( n1ind ) *n1ind = i1;
+  if ( n2ind ) *n2ind = i2;
+  return face;
 }
 
 //=======================================================================
@@ -6588,12 +8161,13 @@ bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode*             theFirst
         vector<const SMDS_MeshNode*> nodes(nbNodes+1);
 
         if(e->IsQuadratic()) {
-          const SMDS_QuadraticFaceOfNodes* F =
-            static_cast<const SMDS_QuadraticFaceOfNodes*>(e);
+          const SMDS_VtkFace* F =
+            dynamic_cast<const SMDS_VtkFace*>(e);
+          if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
           // use special nodes iterator
-          SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
+          SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
           while( anIter->more() ) {
-            nodes[ iNode++ ] = anIter->next();
+            nodes[ iNode++ ] = cast2Node(anIter->next());
           }
         }
         else {
@@ -6763,7 +8337,7 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
     //    links of the free border
     // -------------------------------------------------------------------------
 
-    // 1. Since sewing may brake if there are volumes to split on the side 2,
+    // 1. Since sewing may break if there are volumes to split on the side 2,
     //    we wont move nodes but just compute new coordinates for them
     typedef map<const SMDS_MeshNode*, gp_XYZ> TNodeXYZMap;
     TNodeXYZMap nBordXYZ;
@@ -6861,12 +8435,13 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
         else if ( elem->GetType()==SMDSAbs_Face ) { // --face
           // retrieve all face nodes and find iPrevNode - an index of the prevSideNode
           if(elem->IsQuadratic()) {
-            const SMDS_QuadraticFaceOfNodes* F =
-              static_cast<const SMDS_QuadraticFaceOfNodes*>(elem);
+            const SMDS_VtkFace* F =
+              dynamic_cast<const SMDS_VtkFace*>(elem);
+            if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
             // use special nodes iterator
-            SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
+            SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
             while( anIter->more() ) {
-              nodes[ iNode ] = anIter->next();
+              nodes[ iNode ] = cast2Node(anIter->next());
               if ( nodes[ iNode++ ] == prevSideNode )
                 iPrevNode = iNode - 1;
             }
@@ -7180,12 +8755,13 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement*     theFace,
   vector<const SMDS_MeshNode*> nodes( theFace->NbNodes() );
 
   if(theFace->IsQuadratic()) {
-    const SMDS_QuadraticFaceOfNodes* F =
-      static_cast<const SMDS_QuadraticFaceOfNodes*>(theFace);
+    const SMDS_VtkFace* F =
+      dynamic_cast<const SMDS_VtkFace*>(theFace);
+    if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
     // use special nodes iterator
-    SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
+    SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
     while( anIter->more() ) {
-      const SMDS_MeshNode* n = anIter->next();
+      const SMDS_MeshNode* n = cast2Node(anIter->next());
       if ( n == theBetweenNode1 )
         il1 = iNode;
       else if ( n == theBetweenNode2 )
@@ -7241,12 +8817,13 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement*     theFace,
     bool isFLN = false;
 
     if(theFace->IsQuadratic()) {
-      const SMDS_QuadraticFaceOfNodes* F =
-        static_cast<const SMDS_QuadraticFaceOfNodes*>(theFace);
+      const SMDS_VtkFace* F =
+        dynamic_cast<const SMDS_VtkFace*>(theFace);
+      if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
       // use special nodes iterator
-      SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
+      SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
       while( anIter->more()  && !isFLN ) {
-        const SMDS_MeshNode* n = anIter->next();
+        const SMDS_MeshNode* n = cast2Node(anIter->next());
         poly_nodes[iNode++] = n;
         if (n == nodes[il1]) {
           isFLN = true;
@@ -7259,7 +8836,7 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement*     theFace,
       }
       // add nodes of face starting from last node of link
       while ( anIter->more() ) {
-        poly_nodes[iNode++] = anIter->next();
+        poly_nodes[iNode++] = cast2Node(anIter->next());
       }
     }
     else {
@@ -7302,6 +8879,7 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement*     theFace,
     return;
   }
 
+  SMESHDS_Mesh *aMesh = GetMeshDS();
   if( !theFace->IsQuadratic() ) {
 
     // put aNodesToInsert between theBetweenNode1 and theBetweenNode2
@@ -7352,7 +8930,6 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement*     theFace,
     }
 
     // create new elements
-    SMESHDS_Mesh *aMesh = GetMeshDS();
     int aShapeId = FindShape( theFace );
 
     i1 = 0; i2 = 1;
@@ -7378,8 +8955,16 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement*     theFace,
     newNodes[ 1 ] = linkNodes[ i2 ];
     newNodes[ 2 ] = nodes[ iSplit >= iBestQuad ? i3 : i4 ];
     newNodes[ 3 ] = nodes[ i4 ];
-    aMesh->ChangeElementNodes( theFace, newNodes, iSplit == iBestQuad ? 4 : 3 );
-  } // end if(!theFace->IsQuadratic())
+    //aMesh->ChangeElementNodes( theFace, newNodes, iSplit == iBestQuad ? 4 : 3 );
+    const SMDS_MeshElement* newElem = 0;
+    if (iSplit == iBestQuad)
+      newElem = aMesh->AddFace( newNodes[0], newNodes[1], newNodes[2], newNodes[3] );
+    else
+      newElem = aMesh->AddFace( newNodes[0], newNodes[1], newNodes[2] );
+    myLastCreatedElems.Append(newElem);
+    if ( aShapeId && newElem )
+      aMesh->SetMeshElementOnShape( newElem, aShapeId );
+} // end if(!theFace->IsQuadratic())
   else { // theFace is quadratic
     // we have to split theFace on simple triangles and one simple quadrangle
     int tmp = il1/2;
@@ -7406,7 +8991,6 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement*     theFace,
     //           n4           n6      n5     n4
 
     // create new elements
-    SMESHDS_Mesh *aMesh = GetMeshDS();
     int aShapeId = FindShape( theFace );
 
     int n1,n2,n3;
@@ -7489,9 +9073,9 @@ void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement*     theFace,
       if ( aShapeId && newElem )
         aMesh->SetMeshElementOnShape( newElem, aShapeId );
     }
-    // remove old quadratic face
-    aMesh->RemoveElement(theFace);
   }
+  // remove old face
+  aMesh->RemoveElement(theFace);
 }
 
 //=======================================================================
@@ -7594,7 +9178,7 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
   int nbElem = 0;
   if( !theSm ) return nbElem;
 
-  const bool notFromGroups = false;
+  vector<int> nbNodeInFaces;
   SMDS_ElemIteratorPtr ElemItr = theSm->GetElements();
   while(ElemItr->more())
   {
@@ -7603,16 +9187,14 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
     if( !elem || elem->IsQuadratic() ) continue;
 
     int id = elem->GetID();
+    //MESSAGE("elem " << id);
+    id = 0; // get a free number for new elements
     int nbNodes = elem->NbNodes();
-    vector<const SMDS_MeshNode *> aNds (nbNodes);
-
-    for(int i = 0; i < nbNodes; i++)
-    {
-      aNds[i] = elem->GetNode(i);
-    }
     SMDSAbs_ElementType aType = elem->GetType();
 
-    GetMeshDS()->RemoveFreeElement(elem, theSm, notFromGroups);
+    vector<const SMDS_MeshNode *> nodes (elem->begin_nodes(), elem->end_nodes());
+    if ( elem->GetEntityType() == SMDSEntity_Polyhedra )
+      nbNodeInFaces = static_cast<const SMDS_VtkVolume* >( elem )->GetQuantities();
 
     const SMDS_MeshElement* NewElem = 0;
 
@@ -7620,7 +9202,7 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
     {
     case SMDSAbs_Edge :
       {
-        NewElem = theHelper.AddEdge(aNds[0], aNds[1], id, theForce3d);
+        NewElem = theHelper.AddEdge(nodes[0], nodes[1], id, theForce3d);
         break;
       }
     case SMDSAbs_Face :
@@ -7628,12 +9210,13 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
         switch(nbNodes)
         {
         case 3:
-          NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
+          NewElem = theHelper.AddFace(nodes[0], nodes[1], nodes[2], id, theForce3d);
           break;
         case 4:
-          NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
+          NewElem = theHelper.AddFace(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
           break;
         default:
+          NewElem = theHelper.AddPolygonalFace(nodes, id, theForce3d);
           continue;
         }
         break;
@@ -7643,17 +9226,20 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
         switch(nbNodes)
         {
         case 4:
-          NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
+          NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+          break;
+        case 5:
+          NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], id, theForce3d);
           break;
         case 6:
-          NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], aNds[4], aNds[5], id, theForce3d);
+          NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], nodes[5], id, theForce3d);
           break;
         case 8:
-          NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
-                                        aNds[4], aNds[5], aNds[6], aNds[7], id, theForce3d);
+          NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+                                        nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
           break;
         default:
-          continue;
+          NewElem = theHelper.AddPolyhedralVolume(nodes, nbNodeInFaces, id, theForce3d);
         }
         break;
       }
@@ -7663,7 +9249,11 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
     ReplaceElemInGroups( elem, NewElem, GetMeshDS());
     if( NewElem )
       theSm->AddElement( NewElem );
+
+    GetMeshDS()->RemoveFreeElement(elem, theSm, /*fromGroups=*/false);
   }
+//  if (!GetMeshDS()->isCompacted())
+//    GetMeshDS()->compactMesh();
   return nbElem;
 }
 
@@ -7677,7 +9267,6 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
 
   SMESH_MesherHelper aHelper(*myMesh);
   aHelper.SetIsQuadratic( true );
-  const bool notFromGroups = false;
 
   int nbCheckedElems = 0;
   if ( myMesh->HasShapeToMesh() )
@@ -7705,10 +9294,11 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
       if(edge && !edge->IsQuadratic())
       {
         int id = edge->GetID();
+        //MESSAGE("edge->GetID() " << id);
         const SMDS_MeshNode* n1 = edge->GetNode(0);
         const SMDS_MeshNode* n2 = edge->GetNode(1);
 
-        meshDS->RemoveFreeElement(edge, smDS, notFromGroups);
+        meshDS->RemoveFreeElement(edge, smDS, /*fromGroups=*/false);
 
         const SMDS_MeshEdge* NewEdge = aHelper.AddEdge(n1, n2, id, theForce3d);
         ReplaceElemInGroups( edge, NewEdge, GetMeshDS());
@@ -7722,29 +9312,25 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
 
       int id = face->GetID();
       int nbNodes = face->NbNodes();
-      vector<const SMDS_MeshNode *> aNds (nbNodes);
-
-      for(int i = 0; i < nbNodes; i++)
-      {
-        aNds[i] = face->GetNode(i);
-      }
+      vector<const SMDS_MeshNode *> nodes ( face->begin_nodes(), face->end_nodes());
 
-      meshDS->RemoveFreeElement(face, smDS, notFromGroups);
+      meshDS->RemoveFreeElement(face, smDS, /*fromGroups=*/false);
 
       SMDS_MeshFace * NewFace = 0;
       switch(nbNodes)
       {
       case 3:
-        NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
+        NewFace = aHelper.AddFace(nodes[0], nodes[1], nodes[2], id, theForce3d);
         break;
       case 4:
-        NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
+        NewFace = aHelper.AddFace(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
         break;
       default:
-        continue;
+        NewFace = aHelper.AddPolygonalFace(nodes, id, theForce3d);
       }
       ReplaceElemInGroups( face, NewFace, GetMeshDS());
     }
+    vector<int> nbNodeInFaces;
     SMDS_VolumeIteratorPtr aVolumeItr = meshDS->volumesIterator();
     while(aVolumeItr->more())
     {
@@ -7753,40 +9339,45 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
 
       int id = volume->GetID();
       int nbNodes = volume->NbNodes();
-      vector<const SMDS_MeshNode *> aNds (nbNodes);
+      vector<const SMDS_MeshNode *> nodes (volume->begin_nodes(), volume->end_nodes());
+      if ( volume->GetEntityType() == SMDSEntity_Polyhedra )
+        nbNodeInFaces = static_cast<const SMDS_VtkVolume* >(volume)->GetQuantities();
 
-      for(int i = 0; i < nbNodes; i++)
-      {
-        aNds[i] = volume->GetNode(i);
-      }
-
-      meshDS->RemoveFreeElement(volume, smDS, notFromGroups);
+      meshDS->RemoveFreeElement(volume, smDS, /*fromGroups=*/false);
 
       SMDS_MeshVolume * NewVolume = 0;
       switch(nbNodes)
       {
       case 4:
-        NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
-                                      aNds[3], id, theForce3d );
+        NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
+                                      nodes[3], id, theForce3d );
+        break;
+      case 5:
+        NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
+                                      nodes[3], nodes[4], id, theForce3d);
         break;
       case 6:
-        NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
-                                      aNds[3], aNds[4], aNds[5], id, theForce3d);
+        NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
+                                      nodes[3], nodes[4], nodes[5], id, theForce3d);
         break;
       case 8:
-        NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
-                                      aNds[4], aNds[5], aNds[6], aNds[7], id, theForce3d);
+        NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+                                      nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
         break;
       default:
-        continue;
+        NewVolume = aHelper.AddPolyhedralVolume(nodes, nbNodeInFaces, id, theForce3d);
       }
       ReplaceElemInGroups(volume, NewVolume, meshDS);
     }
   }
-  if ( !theForce3d ) {
-    aHelper.SetSubShape(0); // apply to the whole mesh
+
+  if ( !theForce3d  && !getenv("NO_FixQuadraticElements"))
+  { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
+    aHelper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh
     aHelper.FixQuadraticElements();
   }
+  if (!GetMeshDS()->isCompacted())
+    GetMeshDS()->compactMesh();
 }
 
 //=======================================================================
@@ -7812,8 +9403,8 @@ int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh *    theSm,
     {
       int id = elem->GetID();
       int nbNodes = elem->NbNodes();
-      vector<const SMDS_MeshNode *> aNds, mediumNodes;
-      aNds.reserve( nbNodes );
+      vector<const SMDS_MeshNode *> nodes, mediumNodes;
+      nodes.reserve( nbNodes );
       mediumNodes.reserve( nbNodes );
 
       for(int i = 0; i < nbNodes; i++)
@@ -7823,15 +9414,15 @@ int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh *    theSm,
         if( elem->IsMediumNode( n ) )
           mediumNodes.push_back( n );
         else
-          aNds.push_back( n );
+          nodes.push_back( n );
       }
-      if( aNds.empty() ) continue;
+      if( nodes.empty() ) continue;
       SMDSAbs_ElementType aType = elem->GetType();
 
       //remove old quadratic element
       meshDS->RemoveFreeElement( elem, theSm, notFromGroups );
 
-      SMDS_MeshElement * NewElem = AddElement( aNds, aType, false, id );
+      SMDS_MeshElement * NewElem = AddElement( nodes, aType, false, id );
       ReplaceElemInGroups(elem, NewElem, meshDS);
       if( theSm && NewElem )
         theSm->AddElement( NewElem );
@@ -7841,9 +9432,9 @@ int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh *    theSm,
       for ( ; nIt != mediumNodes.end(); ++nIt ) {
         const SMDS_MeshNode* n = *nIt;
         if ( n->NbInverseElements() == 0 ) {
-          if ( n->GetPosition()->GetShapeId() != theShapeID )
+          if ( n->getshapeId() != theShapeID )
             meshDS->RemoveFreeNode( n, meshDS->MeshElements
-                                    ( n->GetPosition()->GetShapeId() ));
+                                    ( n->getshapeId() ));
           else
             meshDS->RemoveFreeNode( n, theSm );
         }
@@ -7914,12 +9505,13 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
   // 1. Build set of faces representing each side:
   // =======================================================================
   // a. build set of nodes belonging to faces
-  // b. complete set of faces: find missing fices whose nodes are in set of nodes
+  // b. complete set of faces: find missing faces whose nodes are in set of nodes
   // c. create temporary faces representing side of volumes if correspondent
   //    face does not exist
 
   SMESHDS_Mesh* aMesh = GetMeshDS();
-  SMDS_Mesh aTmpFacesMesh;
+  // TODO algoritm not OK with vtkUnstructuredGrid: 2 meshes can't share nodes
+  //SMDS_Mesh aTmpFacesMesh; // try to use the same mesh
   set<const SMDS_MeshElement*> faceSet1, faceSet2;
   set<const SMDS_MeshElement*> volSet1,  volSet2;
   set<const SMDS_MeshNode*>    nodeSet1, nodeSet2;
@@ -7929,6 +9521,7 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
   TIDSortedElemSet * elemSetPtr[] = { &theSide1, &theSide2 };
   int iSide, iFace, iNode;
 
+  list<const SMDS_MeshElement* > tempFaceList;
   for ( iSide = 0; iSide < 2; iSide++ ) {
     set<const SMDS_MeshNode*>    * nodeSet = nodeSetPtr[ iSide ];
     TIDSortedElemSet * elemSet = elemSetPtr[ iSide ];
@@ -7956,7 +9549,7 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
     // -----------------------------------------------------------
     // 1a. Collect nodes of existing faces
     //     and build set of face nodes in order to detect missing
-    //     faces corresponing to sides of volumes
+    //     faces corresponding to sides of volumes
     // -----------------------------------------------------------
 
     set< set <const SMDS_MeshNode*> > setOfFaceNodeSet;
@@ -7980,7 +9573,7 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
         volSet->insert( elem );
     }
     // ------------------------------------------------------------------------------
-    // 1b. Complete set of faces: find missing fices whose nodes are in set of nodes
+    // 1b. Complete set of faces: find missing faces whose nodes are in set of nodes
     // ------------------------------------------------------------------------------
 
     for ( nIt = nodeSet->begin(); nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
@@ -8048,18 +9641,23 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
           if ( !aFreeFace ) {
             // create a temporary face
             if ( nbNodes == 3 ) {
-              aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2] );
+              //aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2] );
+              aFreeFace = aMesh->AddFace( fNodes[0],fNodes[1],fNodes[2] );
             }
             else if ( nbNodes == 4 ) {
-              aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
+              //aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
+              aFreeFace = aMesh->AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
             }
             else {
               vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
-              aFreeFace = aTmpFacesMesh.AddPolygonalFace(poly_nodes);
+              //aFreeFace = aTmpFacesMesh.AddPolygonalFace(poly_nodes);
+              aFreeFace = aMesh->AddPolygonalFace(poly_nodes);
             }
           }
-          if ( aFreeFace )
+          if ( aFreeFace ) {
             freeFaceList.push_back( aFreeFace );
+            tempFaceList.push_back( aFreeFace );
+          }
 
         } // loop on faces of a volume
 
@@ -8089,7 +9687,7 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
               fIt++;
             }
             else
-              freeFaceList.erase( fIt++ ); // here fIt++ occures before erase
+              freeFaceList.erase( fIt++ ); // here fIt++ occurs before erase
           }
           if ( freeFaceList.size() > 1 )
           {
@@ -8173,9 +9771,12 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
 
   if ( faceSet1.size() != faceSet2.size() ) {
     // delete temporary faces: they are in reverseElements of actual nodes
-    SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
-    while ( tmpFaceIt->more() )
-      aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
+//    SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
+//    while ( tmpFaceIt->more() )
+//      aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
+//    list<const SMDS_MeshElement* >::iterator tmpFaceIt = tempFaceList.begin();
+//    for (; tmpFaceIt !=tempFaceList.end(); ++tmpFaceIt)
+//      aMesh->RemoveElement(*tmpFaceIt);
     MESSAGE("Diff nb of faces");
     return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
   }
@@ -8288,10 +9889,11 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
               }
             }
             else { // f->IsQuadratic()
-              const SMDS_QuadraticFaceOfNodes* F =
-                static_cast<const SMDS_QuadraticFaceOfNodes*>(f);
+              const SMDS_VtkFace* F =
+                dynamic_cast<const SMDS_VtkFace*>(f);
+              if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
               // use special nodes iterator
-              SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
+              SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
               while ( anIter->more() ) {
                 const SMDS_MeshNode* n =
                   static_cast<const SMDS_MeshNode*>( anIter->next() );
@@ -8430,9 +10032,12 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
   // ====================================================================
 
   // delete temporary faces: they are in reverseElements of actual nodes
-  SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
-  while ( tmpFaceIt->more() )
-    aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
+//  SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
+//  while ( tmpFaceIt->more() )
+//    aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
+//  list<const SMDS_MeshElement* >::iterator tmpFaceIt = tempFaceList.begin();
+//  for (; tmpFaceIt !=tempFaceList.end(); ++tmpFaceIt)
+//    aMesh->RemoveElement(*tmpFaceIt);
 
   if ( aResult != SEW_OK)
     return aResult;
@@ -8466,7 +10071,21 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
         //         elemIDsToRemove.push_back( e->GetID() );
         //       else
         if ( nbReplaced )
-          aMesh->ChangeElementNodes( e, & nodes[0], nbNodes );
+          {
+            SMDSAbs_ElementType etyp = e->GetType();
+            SMDS_MeshElement* newElem = this->AddElement(nodes, etyp, false);
+            if (newElem)
+              {
+                myLastCreatedElems.Append(newElem);
+                AddToSameGroups(newElem, e, aMesh);
+                int aShapeId = e->getshapeId();
+                if ( aShapeId )
+                  {
+                    aMesh->SetMeshElementOnShape( newElem, aShapeId );
+                  }
+              }
+            aMesh->RemoveElement(e);
+          }
       }
     }
 
@@ -8648,6 +10267,7 @@ SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
   return SEW_OK;
 }
 
+//================================================================================
 /*!
   \brief Creates a hole in a mesh by doubling the nodes of some particular elements
   \param theElems - the list of elements (edges or faces) to be replicated
@@ -8657,6 +10277,8 @@ SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
   replicated nodes should be associated to.
   \return TRUE if operation has been completed successfully, FALSE otherwise
 */
+//================================================================================
+
 bool SMESH_MeshEditor::DoubleNodes( const TIDSortedElemSet& theElems,
                                     const TIDSortedElemSet& theNodesNot,
                                     const TIDSortedElemSet& theAffectedElems )
@@ -8680,6 +10302,7 @@ bool SMESH_MeshEditor::DoubleNodes( const TIDSortedElemSet& theElems,
   return res;
 }
 
+//================================================================================
 /*!
   \brief Creates a hole in a mesh by doubling the nodes of some particular elements
   \param theMeshDS - mesh instance
@@ -8689,6 +10312,8 @@ bool SMESH_MeshEditor::DoubleNodes( const TIDSortedElemSet& theElems,
   \param theIsDoubleElem - flag os to replicate element or modify
   \return TRUE if operation has been completed successfully, FALSE otherwise
 */
+//================================================================================
+
 bool SMESH_MeshEditor::doubleNodes( SMESHDS_Mesh*     theMeshDS,
                                     const TIDSortedElemSet& theElems,
                                     const TIDSortedElemSet& theNodesNot,
@@ -8696,6 +10321,7 @@ bool SMESH_MeshEditor::doubleNodes( SMESHDS_Mesh*     theMeshDS,
                                     const SMDS_MeshNode* >& theNodeNodeMap,
                                     const bool theIsDoubleElem )
 {
+  MESSAGE("doubleNodes");
   // iterate on through element and duplicate them (by nodes duplication)
   bool res = false;
   TIDSortedElemSet::const_iterator elemItr = theElems.begin();
@@ -8724,53 +10350,191 @@ bool SMESH_MeshEditor::doubleNodes( SMESHDS_Mesh*     theMeshDS,
         theNodeNodeMap[ aCurrNode ] = aNewNode;
         myLastCreatedNodes.Append( aNewNode );
       }
-      isDuplicate |= (aCurrNode == aNewNode);
+      isDuplicate |= (aCurrNode != aNewNode);
       newNodes[ ind++ ] = aNewNode;
     }
     if ( !isDuplicate )
       continue;
 
     if ( theIsDoubleElem )
-      myLastCreatedElems.Append( AddElement(newNodes, anElem->GetType(), anElem->IsPoly()) );
+      AddElement(newNodes, anElem->GetType(), anElem->IsPoly());
     else
+      {
+      MESSAGE("ChangeElementNodes");
       theMeshDS->ChangeElementNodes( anElem, &newNodes[ 0 ], anElem->NbNodes() );
-
+      }
     res = true;
   }
   return res;
 }
 
+//================================================================================
 /*!
-  \brief Check if element located inside shape
-  \return TRUE if IN or ON shape, FALSE otherwise
+  \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
+  \return TRUE if operation has been completed successfully, FALSE otherwise
 */
+//================================================================================
 
-static bool isInside(const SMDS_MeshElement* theElem,
-                     BRepClass3d_SolidClassifier& theBsc3d,
-                     const double theTol)
+bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes, 
+                                    const std::list< int >& theListOfModifiedElems )
 {
-  gp_XYZ centerXYZ (0, 0, 0);
-  SMDS_ElemIteratorPtr aNodeItr = theElem->nodesIterator();
-  while (aNodeItr->more())
+  MESSAGE("DoubleNodes");
+  myLastCreatedElems.Clear();
+  myLastCreatedNodes.Clear();
+
+  if ( theListOfNodes.size() == 0 )
+    return false;
+
+  SMESHDS_Mesh* aMeshDS = GetMeshDS();
+  if ( !aMeshDS )
+    return false;
+
+  // iterate through nodes and duplicate them
+
+  std::map< const SMDS_MeshNode*, const SMDS_MeshNode* > anOldNodeToNewNode;
+
+  std::list< int >::const_iterator aNodeIter;
+  for ( aNodeIter = theListOfNodes.begin(); aNodeIter != theListOfNodes.end(); ++aNodeIter )
+  {
+    int aCurr = *aNodeIter;
+    SMDS_MeshNode* aNode = (SMDS_MeshNode*)aMeshDS->FindNode( aCurr );
+    if ( !aNode )
+      continue;
+
+    // duplicate node
+
+    const SMDS_MeshNode* aNewNode = aMeshDS->AddNode( aNode->X(), aNode->Y(), aNode->Z() );
+    if ( aNewNode )
+    {
+      anOldNodeToNewNode[ aNode ] = aNewNode;
+      myLastCreatedNodes.Append( aNewNode );
+    }
+  }
+
+  // Create map of new nodes for modified elements
+
+  std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> > anElemToNodes;
+
+  std::list< int >::const_iterator anElemIter;
+  for ( anElemIter = theListOfModifiedElems.begin(); 
+        anElemIter != theListOfModifiedElems.end(); ++anElemIter )
+  {
+    int aCurr = *anElemIter;
+    SMDS_MeshElement* anElem = (SMDS_MeshElement*)aMeshDS->FindElement( aCurr );
+    if ( !anElem )
+      continue;
+
+    vector<const SMDS_MeshNode*> aNodeArr( anElem->NbNodes() );
+
+    SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
+    int ind = 0;
+    while ( anIter->more() ) 
+    { 
+      SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next();
+      if ( aCurr && anOldNodeToNewNode.find( aCurrNode ) != anOldNodeToNewNode.end() )
+      {
+        const SMDS_MeshNode* aNewNode = anOldNodeToNewNode[ aCurrNode ];
+        aNodeArr[ ind++ ] = aNewNode;
+      }
+      else
+        aNodeArr[ ind++ ] = aCurrNode;
+    }
+    anElemToNodes[ anElem ] = aNodeArr;
+  }
+
+  // Change nodes of elements  
+
+  std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> >::iterator
+    anElemToNodesIter = anElemToNodes.begin();
+  for ( ; anElemToNodesIter != anElemToNodes.end(); ++anElemToNodesIter )
+  {
+    const SMDS_MeshElement* anElem = anElemToNodesIter->first;
+    vector<const SMDS_MeshNode*> aNodeArr = anElemToNodesIter->second;
+    if ( anElem )
+      {
+      MESSAGE("ChangeElementNodes");
+      aMeshDS->ChangeElementNodes( anElem, &aNodeArr[ 0 ], anElem->NbNodes() );
+      }
+  }
+
+  return true;
+}
+
+namespace {
+
+  //================================================================================
+  /*!
+  \brief Check if element located inside shape
+  \return TRUE if IN or ON shape, FALSE otherwise
+  */
+  //================================================================================
+
+  template<class Classifier>
+  bool isInside(const SMDS_MeshElement* theElem,
+                Classifier&             theClassifier,
+                const double            theTol)
   {
-    SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
-    centerXYZ += gp_XYZ(aNode->X(), aNode->Y(), aNode->Z());
+    gp_XYZ centerXYZ (0, 0, 0);
+    SMDS_ElemIteratorPtr aNodeItr = theElem->nodesIterator();
+    while (aNodeItr->more())
+      centerXYZ += SMESH_MeshEditor::TNodeXYZ(cast2Node( aNodeItr->next()));
+
+    gp_Pnt aPnt = centerXYZ / theElem->NbNodes();
+    theClassifier.Perform(aPnt, theTol);
+    TopAbs_State aState = theClassifier.State();
+    return (aState == TopAbs_IN || aState == TopAbs_ON );
   }
-  gp_Pnt aPnt(centerXYZ);
-  theBsc3d.Perform(aPnt, theTol);
-  TopAbs_State aState = theBsc3d.State();
-  return (aState == TopAbs_IN || aState == TopAbs_ON );
+
+  //================================================================================
+  /*!
+   * \brief Classifier of the 3D point on the TopoDS_Face
+   *        with interaface suitable for isInside()
+   */
+  //================================================================================
+
+  struct _FaceClassifier
+  {
+    Extrema_ExtPS       _extremum;
+    BRepAdaptor_Surface _surface;
+    TopAbs_State        _state;
+
+    _FaceClassifier(const TopoDS_Face& face):_extremum(),_surface(face),_state(TopAbs_OUT)
+    {
+      _extremum.Initialize( _surface,
+                            _surface.FirstUParameter(), _surface.LastUParameter(),
+                            _surface.FirstVParameter(), _surface.LastVParameter(),
+                            _surface.Tolerance(), _surface.Tolerance() );
+    }
+    void Perform(const gp_Pnt& aPnt, double theTol)
+    {
+      _state = TopAbs_OUT;
+      _extremum.Perform(aPnt);
+      if ( _extremum.IsDone() )
+        for ( int iSol = 1; iSol <= _extremum.NbExt() && _state == TopAbs_OUT; ++iSol)
+          _state = ( _extremum.Value(iSol) <= theTol ? TopAbs_IN : TopAbs_OUT );
+    }
+    TopAbs_State State() const
+    {
+      return _state;
+    }
+  };
 }
 
+//================================================================================
 /*!
   \brief Creates a hole in a mesh by doubling the nodes of some particular elements
   \param theElems - group of of elements (edges or faces) to be replicated
-  \param theNodesNot - group of nodes not to replicated
+  \param theNodesNot - group of nodes not to replicate
   \param theShape - shape to detect affected elements (element which geometric center
   located on or inside shape).
   The replicated nodes should be associated to affected elements.
   \return TRUE if operation has been completed successfully, FALSE otherwise
 */
+//================================================================================
 
 bool SMESH_MeshEditor::DoubleNodesInRegion( const TIDSortedElemSet& theElems,
                                             const TIDSortedElemSet& theNodesNot,
@@ -8780,8 +10544,17 @@ bool SMESH_MeshEditor::DoubleNodesInRegion( const TIDSortedElemSet& theElems,
     return false;
 
   const double aTol = Precision::Confusion();
-  BRepClass3d_SolidClassifier bsc3d(theShape);
-  bsc3d.PerformInfinitePoint(aTol);
+  auto_ptr< BRepClass3d_SolidClassifier> bsc3d;
+  auto_ptr<_FaceClassifier>              aFaceClassifier;
+  if ( theShape.ShapeType() == TopAbs_SOLID )
+  {
+    bsc3d.reset( new BRepClass3d_SolidClassifier(theShape));;
+    bsc3d->PerformInfinitePoint(aTol);
+  }
+  else if (theShape.ShapeType() == TopAbs_FACE )
+  {
+    aFaceClassifier.reset( new _FaceClassifier(TopoDS::Face(theShape)));
+  }
 
   // iterates on indicated elements and get elements by back references from their nodes
   TIDSortedElemSet anAffected;
@@ -8795,15 +10568,17 @@ bool SMESH_MeshEditor::DoubleNodesInRegion( const TIDSortedElemSet& theElems,
     SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator();
     while ( nodeItr->more() )
     {
-      const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>(nodeItr->next());
+      const SMDS_MeshNode* aNode = cast2Node(nodeItr->next());
       if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() )
         continue;
       SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator();
       while ( backElemItr->more() )
       {
-        SMDS_MeshElement* curElem = (SMDS_MeshElement*)backElemItr->next();
+        const SMDS_MeshElement* curElem = backElemItr->next();
         if ( curElem && theElems.find(curElem) == theElems.end() &&
-             isInside( curElem, bsc3d, aTol ) )
+             ( bsc3d.get() ?
+               isInside( curElem, *bsc3d, aTol ) :
+               isInside( curElem, *aFaceClassifier, aTol )))
           anAffected.insert( curElem );
       }
     }
@@ -8812,10 +10587,200 @@ bool SMESH_MeshEditor::DoubleNodesInRegion( const TIDSortedElemSet& theElems,
 }
 
 /*!
- * \brief Generated skin mesh (containing 2D cells) from 3D mesh
+ * \brief Double nodes on shared faces between groups of volumes and create flat elements on demand.
+ * The list of groups must describe a partition of the mesh volumes.
+ * The nodes of the internal faces at the boundaries of the groups are doubled.
+ * In option, the internal faces are replaced by flat elements.
+ * Triangles are transformed in prisms, and quadrangles in hexahedrons.
+ * @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
+ */
+bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSortedElemSet>& theElems,
+                                                     bool createJointElems)
+{
+  MESSAGE("------------------------------------------------------");
+  MESSAGE("SMESH_MeshEditor::CreateJointElementsOnGroupBoundaries");
+  MESSAGE("------------------------------------------------------");
+
+  SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS();
+  meshDS->BuildDownWardConnectivity(false);
+  CHRONO(50);
+  SMDS_UnstructuredGrid *grid = meshDS->getGrid();
+
+  // --- build the list of faces shared by 2 domains (group of elements), with their domain and volume indexes
+  //     build the list of nodes shared by 2 or more domains, with their domain indexes
+
+  std::map<DownIdType, std::map<int,int>, DownIdCompare> faceDomains; // 2x(id domain --> id volume)
+  std::map<int, std::map<int,int> > nodeDomains; //oldId ->  (domainId -> newId)
+  faceDomains.clear();
+  nodeDomains.clear();
+  std::map<int,int> emptyMap;
+  emptyMap.clear();
+
+  for (int idom = 0; idom < theElems.size(); idom++)
+    {
+
+      // --- build a map (face to duplicate --> volume to modify)
+      //     with all the faces shared by 2 domains (group of elements)
+      //     and corresponding volume of this domain, for each shared face.
+      //     a volume has a face shared by 2 domains if it has a neighbor which is not in is domain.
+
+      const TIDSortedElemSet& domain = theElems[idom];
+      TIDSortedElemSet::const_iterator elemItr = domain.begin();
+      for (; elemItr != domain.end(); ++elemItr)
+        {
+          SMDS_MeshElement* anElem = (SMDS_MeshElement*) *elemItr;
+          if (!anElem)
+            continue;
+          int vtkId = anElem->getVtkId();
+          int neighborsVtkIds[NBMAXNEIGHBORS];
+          int downIds[NBMAXNEIGHBORS];
+          unsigned char downTypes[NBMAXNEIGHBORS];
+          int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId);
+          for (int n = 0; n < nbNeighbors; n++)
+            {
+              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
+                {
+                  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
+                    }
+                }
+            }
+        }
+    }
+
+  MESSAGE("Number of shared faces " << faceDomains.size());
+
+  // --- for each shared face, get the nodes
+  //     for each node, for each domain of the face, create a clone of the node
+
+  std::map<DownIdType, std::map<int,int>, DownIdCompare>::iterator itface = faceDomains.begin();
+  for( ; itface != faceDomains.end();++itface )
+    {
+      DownIdType face = itface->first;
+      std::map<int,int> domvol = itface->second;
+      std::set<int> oldNodes;
+      oldNodes.clear();
+      grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
+      std::set<int>::iterator itn = oldNodes.begin();
+      for (;itn != oldNodes.end(); ++itn)
+        {
+          int oldId = *itn;
+          if (!nodeDomains.count(oldId))
+            nodeDomains[oldId] = emptyMap; // create an empty entry for node
+          std::map<int,int>::iterator itdom = domvol.begin();
+          for(; itdom != domvol.end(); ++itdom)
+            {
+              int idom = itdom->first;
+              if ( nodeDomains[oldId].empty() )
+                nodeDomains[oldId][idom] = oldId; // keep the old node in the first domain
+              else
+                {
+                  double *coords = grid->GetPoint(oldId);
+                  SMDS_MeshNode *newNode = meshDS->AddNode(coords[0], coords[1], coords[2]);
+                  int newId = newNode->getVtkId();
+                  nodeDomains[oldId][idom] = newId; // cloned node for other domains
+                }
+            }
+        }
+    }
+
+  // --- iterate on shared faces (volumes to modify, face to extrude)
+  //     get node id's of the face (id SMDS = id VTK)
+  //     create flat element with old and new nodes if requested
+
+  if (createJointElems)
+    {
+      itface = faceDomains.begin();
+      for( ; itface != faceDomains.end();++itface )
+        {
+          DownIdType face = itface->first;
+          std::set<int> oldNodes;
+          std::set<int>::iterator itn;
+          oldNodes.clear();
+          grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
+          std::map<int,int> localClonedNodeIds;
+
+          std::map<int,int> domvol = itface->second;
+          std::map<int,int>::iterator itdom = domvol.begin();
+          int dom1 = itdom->first;
+          int vtkVolId = itdom->second;
+          itdom++;
+          int dom2 = itdom->first;
+
+          localClonedNodeIds.clear();
+          for (itn = oldNodes.begin(); itn != oldNodes.end(); ++itn)
+            {
+              int oldId = *itn;
+              int refid = oldId;
+              if (nodeDomains[oldId].count(dom1))
+                refid = nodeDomains[oldId][dom1];
+              else
+                MESSAGE("--- problem domain node " << dom1 << " " << oldId);
+              int newid = oldId;
+              if (nodeDomains[oldId].count(dom2))
+                newid = nodeDomains[oldId][dom2];
+              else
+                MESSAGE("--- problem domain node " << dom2 << " " << oldId);
+              localClonedNodeIds[oldId] = newid;
+            }
+          meshDS->extrudeVolumeFromFace(vtkVolId, localClonedNodeIds);
+        }
+    }
+
+  // --- iterate on shared faces (volumes to modify, face to extrude)
+  //     get node id's of the face
+  //     replace old nodes by new nodes in volumes, and update inverse connectivity
+
+  itface = faceDomains.begin();
+  for( ; itface != faceDomains.end();++itface )
+    {
+      DownIdType face = itface->first;
+      std::set<int> oldNodes;
+      std::set<int>::iterator itn;
+      oldNodes.clear();
+      grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
+      std::map<int,int> localClonedNodeIds;
+
+      std::map<int,int> domvol = itface->second;
+      std::map<int,int>::iterator itdom = domvol.begin();
+      for(; itdom != domvol.end(); ++itdom)
+        {
+          int idom = itdom->first;
+          int vtkVolId = itdom->second;
+          localClonedNodeIds.clear();
+          for (itn = oldNodes.begin(); itn != oldNodes.end(); ++itn)
+            {
+              int oldId = *itn;
+              if (nodeDomains[oldId].count(idom))
+                localClonedNodeIds[oldId] = nodeDomains[oldId][idom];
+            }
+          meshDS->ModifyCellNodes(vtkVolId, localClonedNodeIds);
+        }
+    }
+  grid->BuildLinks();
+
+  // TODO replace also old nodes by new nodes in faces and edges
+  CHRONOSTOP(50);
+  counters::stats();
+  return true;
+}
+
+//================================================================================
+/*!
+ * \brief Generates skin mesh (containing 2D cells) from 3D mesh
  * The created 2D mesh elements based on nodes of free faces of boundary volumes
  * \return TRUE if operation has been completed successfully, FALSE otherwise
  */
+//================================================================================
 
 bool SMESH_MeshEditor::Make2DMeshFrom3D()
 {
@@ -8823,46 +10788,225 @@ bool SMESH_MeshEditor::Make2DMeshFrom3D()
   SMESHDS_Mesh* aMesh = GetMeshDS();
   if (!aMesh)
     return false;
-  bool res = false;
+  //bool res = false;
+  int nbFree = 0, nbExisted = 0, nbCreated = 0;
   SMDS_VolumeIteratorPtr vIt = aMesh->volumesIterator();
   while(vIt->more())
   {
     const SMDS_MeshVolume* volume = vIt->next();
     SMDS_VolumeTool vTool( volume );
+    vTool.SetExternalNormal();
     const bool isPoly = volume->IsPoly();
     const bool isQuad = volume->IsQuadratic();
     for ( int iface = 0, n = vTool.NbFaces(); iface < n; iface++ )
     {
       if (!vTool.IsFreeFace(iface))
         continue;
+      nbFree++;
       vector<const SMDS_MeshNode *> nodes;
       int nbFaceNodes = vTool.NbFaceNodes(iface);
       const SMDS_MeshNode** faceNodes = vTool.GetFaceNodes(iface);
-      if (vTool.IsFaceExternal(iface)) 
-      {
-        int inode = 0;
-        for ( ; inode < nbFaceNodes; inode += isQuad ? 2 : 1)
+      int inode = 0;
+      for ( ; inode < nbFaceNodes; inode += isQuad ? 2 : 1)
+        nodes.push_back(faceNodes[inode]);
+      if (isQuad)
+        for ( inode = 1; inode < nbFaceNodes; inode += 2)
           nodes.push_back(faceNodes[inode]);
-        if (isQuad)
-          for ( inode = 1; inode < nbFaceNodes; inode += 2)
-            nodes.push_back(faceNodes[inode]);
+
+      // add new face based on volume nodes
+      if (aMesh->FindFace( nodes ) ) {
+        nbExisted++;
+        continue; // face already exsist
       }
-      else
+      AddElement(nodes, SMDSAbs_Face, isPoly && iface == 1);
+      nbCreated++;
+    }
+  }
+  return ( nbFree==(nbExisted+nbCreated) );
+}
+
+namespace
+{
+  inline const SMDS_MeshNode* getNodeWithSameID(SMESHDS_Mesh* mesh, const SMDS_MeshNode* node)
+  {
+    if ( const SMDS_MeshNode* n = mesh->FindNode( node->GetID() ))
+      return n;
+    return mesh->AddNodeWithID( node->X(),node->Y(),node->Z(), node->GetID() );
+  }
+}
+//================================================================================
+/*!
+ * \brief Creates missing boundary elements
+ *  \param elements - elements whose boundary is to be checked
+ *  \param dimension - defines type of boundary elements to create
+ *  \param group - a group to store created boundary elements in
+ *  \param targetMesh - a mesh to store created boundary elements in
+ *  \param toCopyElements - if true, the checked elements will be copied into the targetMesh
+ *  \param toCopyExistingBondary - if true, not only new but also pre-existing
+ *                                boundary elements will be copied into the targetMesh
+ */
+//================================================================================
+
+void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
+                                        Bnd_Dimension           dimension,
+                                        SMESH_Group*            group/*=0*/,
+                                        SMESH_Mesh*             targetMesh/*=0*/,
+                                        bool                    toCopyElements/*=false*/,
+                                        bool                    toCopyExistingBondary/*=false*/)
+{
+  SMDSAbs_ElementType missType = (dimension == BND_2DFROM3D) ? SMDSAbs_Face : SMDSAbs_Edge;
+  SMDSAbs_ElementType elemType = (dimension == BND_1DFROM2D) ? SMDSAbs_Face : SMDSAbs_Volume;
+  // hope that all elements are of the same type, do not check them all
+  if ( !elements.empty() && (*elements.begin())->GetType() != elemType )
+    throw SALOME_Exception(LOCALIZED("wrong element type"));
+
+  if ( !targetMesh )
+    toCopyElements = toCopyExistingBondary = false;
+
+  SMESH_MeshEditor tgtEditor( targetMesh ? targetMesh : myMesh );
+  SMESHDS_Mesh* aMesh = GetMeshDS(), *tgtMeshDS = tgtEditor.GetMeshDS();
+
+  SMDS_VolumeTool vTool;
+  TIDSortedElemSet emptySet, avoidSet;
+  int inode;
+
+  typedef vector<const SMDS_MeshNode*> TConnectivity;
+
+  SMDS_ElemIteratorPtr eIt;
+  if (elements.empty())
+    eIt = aMesh->elementsIterator(elemType);
+  else
+    eIt = SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() ));
+
+  while (eIt->more())
+  {
+    const SMDS_MeshElement* elem = eIt->next();
+    const int iQuad = elem->IsQuadratic();
+
+    // 1. For an elem, get present bnd elements and connectivities of missing bnd elements
+    vector<const SMDS_MeshElement*> presentBndElems;
+    vector<TConnectivity>           missingBndElems;
+    TConnectivity nodes;
+    if ( vTool.Set(elem) ) // elem is a volume ------------------------------------------
+    {
+      vTool.SetExternalNormal();
+      for ( int iface = 0, n = vTool.NbFaces(); iface < n; iface++ )
       {
-        int inode = nbFaceNodes-1;
-        for ( ; inode >=0; inode -= isQuad ? 2 : 1)
-          nodes.push_back(faceNodes[inode]);
-        if (isQuad)
-          for ( inode = nbFaceNodes-2; inode >=0; inode -= 2)
-            nodes.push_back(faceNodes[inode]);
+        if (!vTool.IsFreeFace(iface))
+          continue;
+        int nbFaceNodes = vTool.NbFaceNodes(iface);
+        const SMDS_MeshNode** nn = vTool.GetFaceNodes(iface);
+        if ( missType == SMDSAbs_Edge ) // boundary edges
+        {
+          nodes.resize( 2+iQuad );
+          for ( int i = 0; i < nbFaceNodes; i += 1+iQuad)
+          {
+            for ( int j = 0; j < nodes.size(); ++j )
+              nodes[j] =nn[i+j];
+            if ( const SMDS_MeshElement* edge =
+                 aMesh->FindElement(nodes,SMDSAbs_Edge,/*noMedium=*/0))
+              presentBndElems.push_back( edge );
+            else
+              missingBndElems.push_back( nodes );
+          }
+        }
+        else // boundary face
+        {
+          nodes.clear();
+          for ( inode = 0; inode < nbFaceNodes; inode += 1+iQuad)
+            nodes.push_back( nn[inode] );
+          if (iQuad)
+            for ( inode = 1; inode < nbFaceNodes; inode += 2)
+              nodes.push_back( nn[inode] );
+
+          if (const SMDS_MeshFace * f = aMesh->FindFace( nodes ) )
+            presentBndElems.push_back( f );
+          else
+            missingBndElems.push_back( nodes );
+        }
+      }
+    }
+    else                     // elem is a face ------------------------------------------
+    {
+      avoidSet.clear(), avoidSet.insert( elem );
+      int nbNodes = elem->NbCornerNodes();
+      nodes.resize( 2 /*+ iQuad*/);
+      for ( int i = 0; i < nbNodes; i++ )
+      {
+        nodes[0] = elem->GetNode(i);
+        nodes[1] = elem->GetNode((i+1)%nbNodes);
+        if ( FindFaceInSet( nodes[0], nodes[1], emptySet, avoidSet))
+          continue; // not free link
+
+        //if ( iQuad )
+        //nodes[2] = elem->GetNode( i + nbNodes );
+        if ( const SMDS_MeshElement* edge =
+             aMesh->FindElement(nodes,SMDSAbs_Edge,/*noMedium=*/true))
+          presentBndElems.push_back( edge );
+        else
+          missingBndElems.push_back( nodes );
       }
+    }
 
-      // add new face based on volume nodes
-      if (aMesh->FindFace( nodes ) )
-        continue; // face already exsist
-      myLastCreatedElems.Append( AddElement(nodes, SMDSAbs_Face, isPoly && iface == 1) );
-      res = true;
+    // 2. Add missing boundary elements
+    if ( targetMesh != myMesh )
+      // instead of making a map of nodes in this mesh and targetMesh,
+      // we create nodes with same IDs. We can renumber them later, if needed
+      for ( int i = 0; i < missingBndElems.size(); ++i )
+      {
+        TConnectivity& srcNodes = missingBndElems[i];
+        TConnectivity  nodes( srcNodes.size() );
+        for ( inode = 0; inode < nodes.size(); ++inode )
+          nodes[inode] = getNodeWithSameID( tgtMeshDS, srcNodes[inode] );
+        tgtEditor.AddElement(nodes, missType, elem->IsPoly() && nodes.size()/(iQuad+1)>4);
+      }
+    else
+      for ( int i = 0; i < missingBndElems.size(); ++i )
+      {
+        TConnectivity&  nodes = missingBndElems[i];
+        tgtEditor.AddElement(nodes, missType, elem->IsPoly() && nodes.size()/(iQuad+1)>4);
+      }
+
+    // 3. Copy present boundary elements
+    if ( toCopyExistingBondary )
+      for ( int i = 0 ; i < presentBndElems.size(); ++i )
+      {
+        const SMDS_MeshElement* e = presentBndElems[i];
+        TConnectivity nodes( e->NbNodes() );
+        for ( inode = 0; inode < nodes.size(); ++inode )
+          nodes[inode] = getNodeWithSameID( tgtMeshDS, e->GetNode(inode) );
+        tgtEditor.AddElement(nodes, missType, e->IsPoly());
+        // leave only missing elements in tgtEditor.myLastCreatedElems
+        tgtEditor.myLastCreatedElems.Remove( tgtEditor.myLastCreatedElems.Size() );
+      }
+  } // loop on given elements
+
+  // 4. Fill group with missing boundary elements
+  if ( group )
+  {
+    if ( SMESHDS_Group* g = dynamic_cast<SMESHDS_Group*>( group->GetGroupDS() ))
+      for ( int i = 0; i < tgtEditor.myLastCreatedElems.Size(); ++i )
+        g->SMDSGroup().Add( tgtEditor.myLastCreatedElems( i+1 ));
+  }
+  tgtEditor.myLastCreatedElems.Clear();
+
+  // 5. Copy given elements
+  if ( toCopyElements )
+  {
+    if (elements.empty())
+      eIt = aMesh->elementsIterator(elemType);
+    else
+      eIt = SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() ));
+    while (eIt->more())
+    {
+      const SMDS_MeshElement* elem = eIt->next();
+      TConnectivity nodes( elem->NbNodes() );
+      for ( inode = 0; inode < nodes.size(); ++inode )
+        nodes[inode] = getNodeWithSameID( tgtMeshDS, elem->GetNode(inode) );
+      tgtEditor.AddElement(nodes, elemType, elem->IsPoly());
+
+      tgtEditor.myLastCreatedElems.Clear();
     }
   }
-  return res;
+  return;
 }