Salome HOME
0020889: EDF 1433 SMESH: SplitHexaToTetra: add the 24 tetras mode
[modules/smesh.git] / src / SMESH / SMESH_MeshEditor.cxx
index 3f496e87fd06bf75dab5081b6b5c711c83d0df5b..b802540f09405bf6c1902001f5da04a41d31a499 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
@@ -19,6 +19,7 @@
 //
 //  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
@@ -124,6 +125,11 @@ SMESH_MeshEditor::AddElement(const vector<const SMDS_MeshNode*> & node,
   int nbnode = node.size();
   SMESHDS_Mesh* mesh = GetMeshDS();
   switch ( type ) {
+  case SMDSAbs_0DElement:
+    if ( nbnode == 1 )
+      if ( ID ) e = mesh->Add0DElementWithID(node[0], ID);
+      else      e = mesh->Add0DElement      (node[0] );
+    break;
   case SMDSAbs_Edge:
     if ( nbnode == 2 )
       if ( ID ) e = mesh->AddEdgeWithID(node[0], node[1], ID);
@@ -245,8 +251,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();
@@ -254,6 +260,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;
@@ -290,6 +297,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
@@ -303,7 +311,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;
 }
 
 //=======================================================================
@@ -1156,159 +1164,537 @@ namespace
 {
   // Methods of splitting volumes into tetra
 
-  const int theHexTo5[5*4] =
+  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] =
     {
-      0, 1, 5, 2,
-      0, 4, 5, 7,
-      0, 3, 7, 2,
-      5, 6, 7, 2,
-      0, 2, 5, 7
+      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[6*4] =
+  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, 5, 2,
-      0, 4, 5, 7,
-      0, 3, 7, 2,
-      5, 6, 7, 2,
-      0, 2, 5, 7
+      0, 1, 2, 4,    0, 2, 3, 4,   -1
     };
-  const int thePyraTo2[2*4] =
+  const int thePyraTo2_2[2*4+1] =
     {
-      0, 1, 2, 4,
-      0, 2, 3, 4
+      1, 2, 3, 4,    1, 3, 0, 4,   -1
     };
+  const int* thePyraTo2[2] = { thePyraTo2_1, thePyraTo2_2 };
 
-  const int thePentaTo8[8*4] =
+  const int thePentaTo3_1[3*4+1] =
     {
-      0, 1, 2, 6,
-      3, 5, 4, 6,
-      0, 3, 4, 6,
-      0, 4, 1, 6,
-      1, 4, 5, 6,
-      1, 5, 2, 6,
-      2, 5, 3, 6,
-      2, 3, 0, 6
+      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;
-    bool       _addNode; // additional node is to be created
+    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), _addNode(addNode) {}
+      : _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( const SMDS_MeshElement* vol, const int theMethodFlags)
+  //=======================================================================
+
+  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->GetType() == SMDSAbs_Volume && !vol->IsPoly())
-      switch ( vol->NbNodes() )
+    if ( !vol.Element()->IsPoly() && !is24TetMode )
+    {
+      int nbVariants = 2, nbTet = 0;
+      const int** connVariants = 0;
+      switch ( vol.Element()->GetEntityType() )
       {
-      case 8:
-      case 20:
-        if ( theMethodFlags & SMESH_MeshEditor::HEXA_TO_5 )
-          method = TSplitMethod( 5, theHexTo5 );
+      case SMDSEntity_Hexa:
+      case SMDSEntity_Quad_Hexa:
+        if ( theMethodFlags == SMESH_MeshEditor::HEXA_TO_5 )
+          connVariants = theHexTo5, nbTet = 5;
         else
-          method = TSplitMethod( 6, theHexTo6 );
+          connVariants = theHexTo6, nbTet = 6, nbVariants = 4;
         break;
-      case 5:
-      case 13:
-        method = TSplitMethod( 2, thePyraTo2 );
+      case SMDSEntity_Pyramid:
+      case SMDSEntity_Quad_Pyramid:
+        connVariants = thePyraTo2;  nbTet = 2;
         break;
-      case 6:
-      case 15:
-        method = TSplitMethod( 8, thePentaTo8, /*addNode=*/true );
+      case SMDSEntity_Penta:
+      case SMDSEntity_Quad_Penta:
+        connVariants = thePentaTo3; nbTet = 3; nbVariants = 6;
         break;
-      default:;
+      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)
-// {
-//   // sdt-like iterator on coordinates of nodes of mesh element
-//   typedef SMDS_StdIterator< TNodeXYZ, SMDS_ElemIteratorPtr > NXyzIterator;
-//   NXyzIterator xyzEnd;
-
-//   SMESH_MesherHelper helper( *GetMesh());
-
-//   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 edge
-
-//     if ( (*elem)->IsQuadratic() )
-//     {
-//       // add quadratic links to the helper
-//       SMDS_VolumeTool vol( *elem );
-//       for ( int iF = 0; iF < vol.NbFaces(); ++iF )
-//       {
-//         const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iF );
-//         for ( int iN = 0; iN < vol.NbFaceNodes( iF ); iN += 2)
-//           helper.AddTLinkNode( fNodes[iF], fNodes[iF+2], fNodes[iF+1] );
-//       }
-//       helper.SetIsQuadratic( true );
-//     }
-//     else
-//     {
-//       helper.SetIsQuadratic( false );
-//     }
-
-//     vector<const SMDS_MeshElement* > tetras; // splits of a volume
-
-//     if ( geomType == SMDSEntity_Polyhedra )
-//     {
-//       // Each face of a polyhedron is split into triangles and
-//       // each of triangles and a cell barycenter form a tetrahedron.
-
-//       SMDS_VolumeTool vol( *elem );
-
-//       // make a node at barycenter
-//       gp_XYZ gc = std::accumulate( NXyzIterator((*elem)->nodesIterator()), xyzEnd,gp_XYZ(0,0,0));
-//       gc /= vol.NbNodes();
-//       SMDS_MeshNode* gcNode = GetMeshDS()->AddNode( gc.X(), gc.Y(), gc.Z() );
-
-//       for ( int iF = 0; iF < vol.NbFaces(); ++iF )
-//       {
-//         const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iF );
-//         int nbFNodes = vol.NbFaceNodes( iF );
-//         int nbTria = nbFNodes - 2;
-//         bool extFace =  vol.IsFaceExternal( iF );
-//         SMDS_MeshElement* tet;
-//         for ( int i = 0; i < nbTria; ++i )
-//         {
-//           if ( extFace )
-//             tet = helper.AddVolume( fNodes[0], fNodes[i+1], fNodes[i+2], gcNode );
-//           else
-//             tet = helper.AddVolume( fNodes[0], fNodes[i+2], fNodes[i+1], gcNode );
-//           tetras.push_back( tet );
-//         }
-//       }
-
-//     }
-//     else
-//     {
-
-//       TSplitMethod splitMethod = getSplitMethod( *elem, theMethodFlags );
-//       if ( splitMethod._nbTetra < 1 ) continue;
-
-//       vector<const SMDS_MeshNode*> volNodes( (*elem)->begin_nodes(), (*elem)->end_nodes());
-//     }
-//   }
-// }
+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.back() ) continue;
+          if ( fSubMesh )
+            fSubMesh->AddElement( triangles.back());
+          newElems.Append( triangles.back() );
+        }
+        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
@@ -1351,10 +1737,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,
@@ -1371,6 +1758,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.
@@ -4787,10 +5197,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,
@@ -4818,6 +5235,7 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
     groupPostfix = "translated";
     break;
   case gp_Scale:
+  case gp_CompoundTrsf: // different scale by axis
     groupPostfix = "scaled";
     break;
   default:
@@ -4840,7 +5258,36 @@ 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
+  list<SMDS_MeshNode>          orphanCopy; // copies of orphan nodes
+  vector<const SMDS_MeshNode*> orphanNode; // original orphan nodes
+
+  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_MeshElementIDFactory idFactory;
+    SMDS_NodeIteratorPtr nIt = aMesh->nodesIterator();
+    while ( nIt->more() )
+    {
+      const SMDS_MeshNode* node = nIt->next();
+      if ( node->NbInverseElements() == 0 && !theElems.insert( node ).second )
+      {
+        // node was not inserted into theElems because an element with the same ID
+        // is already there. As a work around we insert a copy of node with
+        // an ID = -<index in orphanNode>
+        orphanCopy.push_back( *node ); // copy node
+        SMDS_MeshNode* nodeCopy = &orphanCopy.back();
+        int uniqueID = -orphanNode.size();
+        orphanNode.push_back( node );
+        idFactory.BindID( uniqueID, nodeCopy );
+        theElems.insert( nodeCopy );
+      }
+    }
+  }
+  // loop on theElems to transorm nodes
   TIDSortedElemSet::iterator itElem;
   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
     const SMDS_MeshElement* elem = *itElem;
@@ -4851,8 +5298,10 @@ 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() );
+      if ( node->GetID() < 0 )
+        node = orphanNode[ -node->GetID() ];
+      // check if a node has been already transformed
       pair<TNodeNodeMap::iterator,bool> n2n_isnew =
         nodeMap.insert( make_pair ( node, node ));
       if ( !n2n_isnew.second )
@@ -5103,331 +5552,6 @@ 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)
-{
-  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_PolyhedralVolumeOfNodes* aPolyedre =
-            dynamic_cast<const SMDS_PolyhedralVolumeOfNodes*>( 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
@@ -5562,24 +5686,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);
 }
 
 
@@ -5599,9 +5720,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() );
     }
@@ -5653,12 +5774,13 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
       treeList.push_back( myOctreeNode );
 
       SMDS_MeshNode 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() );
@@ -5949,6 +6071,9 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
                                   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);
@@ -5957,7 +6082,7 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
   {
     return _outerFaces.empty() || _outerFaces.count(face);
   }
-  struct TInters //!< data of intersection of the line and the mesh face used in GetPointState()
+  struct TInters //!< data of intersection of the line and the mesh face (used in GetPointState())
   {
     const SMDS_MeshElement* _face;
     gp_Vec                  _faceNorm;
@@ -6458,6 +6583,26 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
   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 );
+  }
+  TIDSortedElemSet suspectFaces; // elements possibly intersecting the line
+  _ebbTree->getElementsNearLine( line, suspectFaces );
+  foundElems.assign( suspectFaces.begin(), suspectFaces.end());
+}
+
 //=======================================================================
 /*!
  * \brief Return SMESH_ElementSearcher
@@ -8541,7 +8686,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())
   {
@@ -8551,15 +8696,13 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
 
     int id = elem->GetID();
     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_PolyhedralVolumeOfNodes* >( elem )->GetQuanities();
+
+    GetMeshDS()->RemoveFreeElement(elem, theSm, /*fromGroups=*/false);
 
     const SMDS_MeshElement* NewElem = 0;
 
@@ -8567,7 +8710,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 :
@@ -8575,12 +8718,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;
@@ -8590,20 +8734,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(aNds[0], aNds[1], aNds[2], aNds[3], aNds[4], id, theForce3d);
+          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;
       }
@@ -8627,7 +8771,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() )
@@ -8658,7 +8801,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
         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());
@@ -8672,29 +8815,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())
     {
@@ -8703,42 +8842,41 @@ 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_PolyhedralVolumeOfNodes* >(volume)->GetQuanities();
 
-      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(aNds[0], aNds[1], aNds[2],
-                                      aNds[3], aNds[4], id, theForce3d);
+        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();
   }
 }
@@ -8766,8 +8904,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++)
@@ -8777,15 +8915,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 );