Salome HOME
BUG: EDF 2655: Low performance of hexa to tetra splitting
[modules/smesh.git] / src / SMESH / SMESH_MeshEditor.cxx
index c62006a0eb2508d7ec9df9ba2f59518ba034dc4c..1d5698d249bc67e59ed207f38ab86ebc735811ca 100644 (file)
 using namespace std;
 using namespace SMESH::Controls;
 
-typedef SMDS_SetIterator< SMDS_pElement, TIDSortedElemSet::const_iterator> TSetIterator;
+namespace
+{
+  template < class ELEM_SET >
+  SMDS_ElemIteratorPtr elemSetIterator( const ELEM_SET& elements )
+  {
+    typedef SMDS_SetIterator
+      < SMDS_pElement, typename ELEM_SET::const_iterator> TSetIterator;
+    return SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() ));
+  }
+}
 
 //=======================================================================
 //function : SMESH_MeshEditor
@@ -416,12 +425,21 @@ int SMESH_MeshEditor::Remove (const list< int >& theIDs,
 void SMESH_MeshEditor::Create0DElementsOnAllNodes( const TIDSortedElemSet& elements,
                                                    TIDSortedElemSet&       all0DElems )
 {
-  typedef SMDS_SetIterator<const SMDS_MeshElement*, TIDSortedElemSet::const_iterator> TSetIterator;
   SMDS_ElemIteratorPtr elemIt;
+  vector< const SMDS_MeshElement* > allNodes;
   if ( elements.empty() )
+  {
+    allNodes.reserve( GetMeshDS()->NbNodes() );
     elemIt = GetMeshDS()->elementsIterator( SMDSAbs_Node );
+    while ( elemIt->more() )
+      allNodes.push_back( elemIt->next() );
+
+    elemIt = elemSetIterator( allNodes );
+  }
   else
-    elemIt = SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() ));
+  {
+    elemIt = elemSetIterator( elements );
+  }
 
   while ( elemIt->more() )
   {
@@ -1254,8 +1272,6 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet &                   theElems,
   myLastCreatedElems.Clear();
   myLastCreatedNodes.Clear();
 
-  MESSAGE( "::QuadToTri()" );
-
   if ( !theCrit.get() )
     return false;
 
@@ -1347,6 +1363,130 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet &                   theElems,
   return true;
 }
 
+//=======================================================================
+/*!
+ * \brief Split each of given quadrangles into 4 triangles.
+ * \param theElems - The faces to be splitted. If empty all faces are split.
+ */
+//=======================================================================
+
+void SMESH_MeshEditor::QuadTo4Tri (TIDSortedElemSet & theElems)
+{
+  myLastCreatedElems.Clear();
+  myLastCreatedNodes.Clear();
+
+  SMESH_MesherHelper helper( *GetMesh() );
+  helper.SetElementsOnShape( true );
+
+  SMDS_ElemIteratorPtr faceIt;
+  if ( theElems.empty() ) faceIt = GetMeshDS()->elementsIterator(SMDSAbs_Face);
+  else                    faceIt = elemSetIterator( theElems );
+
+  bool   checkUV;
+  gp_XY  uv [9]; uv[8] = gp_XY(0,0);
+  gp_XYZ xyz[9];
+  vector< const SMDS_MeshNode* > nodes;
+  SMESHDS_SubMesh*               subMeshDS;
+  TopoDS_Face                    F;
+  Handle(Geom_Surface)           surface;
+  TopLoc_Location                loc;
+
+  while ( faceIt->more() )
+  {
+    const SMDS_MeshElement* quad = faceIt->next();
+    if ( !quad || quad->NbCornerNodes() != 4 )
+      continue;
+
+    // get a surface the quad is on
+
+    if ( quad->getshapeId() < 1 )
+    {
+      F.Nullify();
+      helper.SetSubShape( 0 );
+      subMeshDS = 0;
+    }
+    else if ( quad->getshapeId() != helper.GetSubShapeID() )
+    {
+      helper.SetSubShape( quad->getshapeId() );
+      if ( !helper.GetSubShape().IsNull() &&
+           helper.GetSubShape().ShapeType() == TopAbs_FACE )
+      {
+        F = TopoDS::Face( helper.GetSubShape() );
+        surface = BRep_Tool::Surface( F, loc );
+        subMeshDS = GetMeshDS()->MeshElements( quad->getshapeId() );
+      }
+      else
+      {
+        helper.SetSubShape( 0 );
+        subMeshDS = 0;
+      }
+    }
+
+    // create a central node
+
+    const SMDS_MeshNode* nCentral;
+    nodes.assign( quad->begin_nodes(), quad->end_nodes() );
+
+    if ( nodes.size() == 9 )
+    {
+      nCentral = nodes.back();
+    }
+    else
+    {
+      size_t iN = 0;
+      if ( F.IsNull() )
+      {
+        for ( ; iN < nodes.size(); ++iN )
+          xyz[ iN ] = SMESH_TNodeXYZ( nodes[ iN ] );
+
+        for ( ; iN < 8; ++iN ) // mid-side points of a linear qudrangle
+          xyz[ iN ] = 0.5 * ( xyz[ iN - 4 ] + xyz[( iN - 3 )%4 ] );
+
+        xyz[ 8 ] = helper.calcTFI( 0.5, 0.5,
+                                   xyz[0], xyz[1], xyz[2], xyz[3],
+                                   xyz[4], xyz[5], xyz[6], xyz[7] );
+      }
+      else
+      {
+        for ( ; iN < nodes.size(); ++iN )
+          uv[ iN ] = helper.GetNodeUV( F, nodes[iN], nodes[(iN+2)%4], &checkUV );
+
+        for ( ; iN < 8; ++iN ) // UV of mid-side points of a linear qudrangle
+          uv[ iN ] = helper.GetMiddleUV( surface, uv[ iN - 4 ], uv[( iN - 3 )%4 ] );
+
+        uv[ 8 ] = helper.calcTFI( 0.5, 0.5,
+                                  uv[0], uv[1], uv[2], uv[3],
+                                  uv[4], uv[5], uv[6], uv[7] );
+
+        gp_Pnt p = surface->Value( uv[8].X(), uv[8].Y() ).Transformed( loc );
+        xyz[ 8 ] = p.XYZ();
+      }
+
+      nCentral = helper.AddNode( xyz[8].X(), xyz[8].Y(), xyz[8].Z(), /*id=*/0,
+                                 uv[8].X(), uv[8].Y() );
+      myLastCreatedNodes.Append( nCentral );
+    }
+
+    // create 4 triangles
+
+    GetMeshDS()->RemoveFreeElement( quad, subMeshDS, /*fromGroups=*/false );
+    
+    helper.SetIsQuadratic  ( nodes.size() > 4 );
+    helper.SetIsBiQuadratic( nodes.size() == 9 );
+    if ( helper.GetIsQuadratic() )
+      helper.AddTLinks( static_cast< const SMDS_MeshFace*>( quad ));
+
+    for ( int i = 0; i < 4; ++i )
+    {
+      SMDS_MeshElement* tria = helper.AddFace( nodes[ i ],
+                                               nodes[(i+1)%4],
+                                               nCentral );
+      ReplaceElemInGroups( tria, quad, GetMeshDS() );
+      myLastCreatedElems.Append( tria );
+    }
+  }
+}
+
 //=======================================================================
 //function : BestSplit
 //purpose  : Find better diagonal for cutting.
@@ -1759,6 +1899,20 @@ namespace
   };
 } // namespace
 
+class TElemToDelete
+{
+public:
+  TElemToDelete(const SMDS_MeshElement* theElem, SMESHDS_SubMesh* theSubMesh)
+  {
+    elem = theElem;
+    subMesh = theSubMesh;
+  }
+  const SMDS_MeshElement* Elem() const {return elem;}
+  SMESHDS_SubMesh* Submesh() {return subMesh;}
+  const SMDS_MeshElement* elem;
+  SMESHDS_SubMesh* subMesh;
+};
+
 //=======================================================================
 //function : SplitVolumesIntoTetra
 //purpose  : Split volume elements into tetrahedra.
@@ -1784,6 +1938,7 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
   double bc[3];
 
   TIDSortedElemSet::const_iterator elem = theElems.begin();
+  std::vector<TElemToDelete> elem_to_delete;
   for ( ; elem != theElems.end(); ++elem )
   {
     if ( (*elem)->GetType() != SMDSAbs_Volume )
@@ -1947,11 +2102,26 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
         }
         ReplaceElemInGroups( face, triangles, GetMeshDS() );
         GetMeshDS()->RemoveFreeElement( face, fSubMesh, /*fromGroups=*/false );
+//         TElemToDelete faceToDelete(face, fSubMesh);
+//         elem_to_delete.push_back(faceToDelete); 
       }
 
     } // loop on volume faces to split them into triangles
 
-    GetMeshDS()->RemoveFreeElement( *elem, subMesh, /*fromGroups=*/false );
+//     GetMeshDS()->RemoveFreeElement( *elem, subMesh, /*fromGroups=*/false );
+
+    // rnc : don't delete the elem here because it results in a mesh with a free
+    // ID at the beginning of the ID list. The first tetra is then inserted in O(1)
+    // but the second one is inserted in O(n), then the whole procedure has almost a O(n^2)
+    // complexity. If all elements to remove are stored and removed after tetra creation 
+    // we get a O(n) complexity for the whole procedure. 
+    // The memory cost is at worst a 6*n*constant memory occupation (where n is the number of elements) 
+    // before deletion of the hexas and then 5*n*constant instead of a maximum of 5*n*constant.
+    // So there is a transient 1/5*(memory occupation) additional cost.
+
+    // Store the elements to delete
+    TElemToDelete elemToDelete(*elem, subMesh);
+    elem_to_delete.push_back(elemToDelete);
 
     if ( geomType == SMDSEntity_TriQuad_Hexa )
     {
@@ -1962,6 +2132,13 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
     }
   } // loop on volumes to split
 
+  // Delete stored elements
+  std::vector<TElemToDelete>::iterator it;
+  for( it = elem_to_delete.begin(); it!= elem_to_delete.end(); it++)
+  {
+    GetMeshDS()->RemoveFreeElement( it->Elem(), it->Submesh(), /*fromGroups=*/false );
+  }
+  
   myLastCreatedNodes = newNodes;
   myLastCreatedElems = newElems;
 }
@@ -5004,7 +5181,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
 
   const TopoDS_Shape& aS = theTrack->GetShapeToMesh();
 
-  if( aS == SMESH_Mesh::PseudoShape() ) {
+  if ( !theTrack->HasShapeToMesh() ) {
     //Mesh without shape
     const SMDS_MeshNode* currentNode = NULL;
     const SMDS_MeshNode* prevNode = theN1;
@@ -5127,16 +5304,8 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
     if ( BRep_Tool::Degenerated( aTrackEdge ) )
       return EXTR_BAD_PATH_SHAPE;
     TopExp::Vertices( aTrackEdge, aV1, aV2 );
-    const SMDS_MeshNode* aN1 = 0;
-    const SMDS_MeshNode* aN2 = 0;
-    if ( theTrack->GetSubMesh( aV1 ) && theTrack->GetSubMesh( aV1 )->GetSubMeshDS() ) {
-      aItN = theTrack->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
-      aN1 = aItN->next();
-    }
-    if ( theTrack->GetSubMesh( aV2 ) && theTrack->GetSubMesh( aV2 )->GetSubMeshDS() ) {
-      aItN = theTrack->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
-      aN2 = aItN->next();
-    }
+    const SMDS_MeshNode* aN1 = SMESH_Algo::VertexNode( aV1, pMeshDS );
+    const SMDS_MeshNode* aN2 = SMESH_Algo::VertexNode( aV2, pMeshDS );
     // starting node must be aN1 or aN2
     if ( !( aN1 == theN1 || aN2 == theN1 ) )
       return EXTR_BAD_STARTING_NODE;
@@ -5184,17 +5353,8 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
         if ( aVprev.IsNull() ) {
           // if previous vertex is not yet defined, it means that we in the beginning of wire
           // and we have to find initial vertex corresponding to starting node theN1
-          const SMDS_MeshNode* aN1 = 0;
-          const SMDS_MeshNode* aN2 = 0;
-
-          if ( locTrack->GetFather()->GetSubMesh(aV1) && locTrack->GetFather()->GetSubMesh(aV1)->GetSubMeshDS() ) {
-            aItN = locTrack->GetFather()->GetSubMesh(aV1)->GetSubMeshDS()->GetNodes();
-            aN1 = aItN->next();
-          }
-          if ( locTrack->GetFather()->GetSubMesh(aV2) && locTrack->GetFather()->GetSubMesh(aV2)->GetSubMeshDS() ) {
-            aItN = locTrack->GetFather()->GetSubMesh(aV2)->GetSubMeshDS()->GetNodes();
-            aN2 = aItN->next();
-          }
+          const SMDS_MeshNode* aN1 = SMESH_Algo::VertexNode( aV1, pMeshDS );
+          const SMDS_MeshNode* aN2 = SMESH_Algo::VertexNode( aV2, pMeshDS );
           // starting node must be aN1 or aN2
           aN1isOK = ( aN1 && aN1 == theN1 );
           aN2isOK = ( aN2 && aN2 == theN1 );
@@ -5228,27 +5388,21 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
       }
     }
     list< list<SMESH_MeshEditor_PathPoint> >::iterator itLLPP = LLPPs.begin();
-    list<SMESH_MeshEditor_PathPoint> firstList = *itLLPP;
-    list<SMESH_MeshEditor_PathPoint>::iterator itPP = firstList.begin();
-    for(; itPP!=firstList.end(); itPP++) {
-      fullList.push_back( *itPP );
-    }
+    list<SMESH_MeshEditor_PathPoint>& firstList = *itLLPP;
+    fullList.splice( fullList.end(), firstList );
+
     SMESH_MeshEditor_PathPoint PP1 = fullList.back();
     fullList.pop_back();
     itLLPP++;
     for(; itLLPP!=LLPPs.end(); itLLPP++) {
-      list<SMESH_MeshEditor_PathPoint> currList = *itLLPP;
-      itPP = currList.begin();
+      list<SMESH_MeshEditor_PathPoint>& currList = *itLLPP;
       SMESH_MeshEditor_PathPoint PP2 = currList.front();
       gp_Dir D1 = PP1.Tangent();
       gp_Dir D2 = PP2.Tangent();
       gp_Dir Dnew( ( D1.XYZ() + D2.XYZ() ) / 2 );
       PP1.SetTangent(Dnew);
       fullList.push_back(PP1);
-      itPP++;
-      for(; itPP!=currList.end(); itPP++) {
-        fullList.push_back( *itPP );
-      }
+      fullList.splice( fullList.end(), currList, ++currList.begin(), currList.end() );
       PP1 = fullList.back();
       fullList.pop_back();
     }
@@ -5270,9 +5424,9 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
 //purpose  : auxilary for ExtrusionAlongTrack
 //=======================================================================
 SMESH_MeshEditor::Extrusion_Error
-SMESH_MeshEditor::MakeEdgePathPoints(std::list<double>& aPrms,
-                                     const TopoDS_Edge& aTrackEdge,
-                                     bool FirstIsStart,
+SMESH_MeshEditor::MakeEdgePathPoints(std::list<double>&                aPrms,
+                                     const TopoDS_Edge&                aTrackEdge,
+                                     bool                              FirstIsStart,
                                      list<SMESH_MeshEditor_PathPoint>& LPP)
 {
   Standard_Real aTx1, aTx2, aL2, aTolVec, aTolVec2;
@@ -5325,63 +5479,43 @@ SMESH_MeshEditor::MakeEdgePathPoints(std::list<double>& aPrms,
 //purpose  : auxilary for ExtrusionAlongTrack
 //=======================================================================
 SMESH_MeshEditor::Extrusion_Error
-SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet&  theElements,
+SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet&                 theElements,
                                    list<SMESH_MeshEditor_PathPoint>& fullList,
-                                   const bool theHasAngles,
-                                   list<double>& theAngles,
-                                   const bool theLinearVariation,
-                                   const bool theHasRefPoint,
-                                   const gp_Pnt& theRefPoint,
-                                   const bool theMakeGroups)
+                                   const bool                        theHasAngles,
+                                   list<double>&                     theAngles,
+                                   const bool                        theLinearVariation,
+                                   const bool                        theHasRefPoint,
+                                   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);
+  const int aNbTP = fullList.size();
   // Angles
-  if( theHasAngles && theAngles.size()>0 && theLinearVariation ) {
+  if( theHasAngles && !theAngles.empty() && theLinearVariation )
     LinearAngleVariation(aNbTP-1, theAngles);
-  }
-  vector<double> aAngles( aNbTP );
-  int j = 0;
-  for(; j<aNbTP; ++j) {
-    aAngles[j] = 0.;
-  }
-  if ( theHasAngles ) {
-    double anAngle;;
-    std::list<double>::iterator aItD = theAngles.begin();
-    for ( j=1; (aItD != theAngles.end()) && (j<aNbTP); ++aItD, ++j ) {
-      anAngle = *aItD;
-      aAngles[j] = anAngle;
-    }
-  }
   // fill vector of path points with angles
-  //aPPs.resize(fullList.size());
-  j = -1;
+  vector<SMESH_MeshEditor_PathPoint> aPPs;
   list<SMESH_MeshEditor_PathPoint>::iterator itPP = fullList.begin();
-  for(; itPP!=fullList.end(); itPP++) {
-    j++;
-    SMESH_MeshEditor_PathPoint PP = *itPP;
-    PP.SetAngle(aAngles[j]);
-    aPPs[j] = PP;
+  list<double>::iterator                 itAngles = theAngles.begin();
+  aPPs.push_back( *itPP++ );
+  for( ; itPP != fullList.end(); itPP++) {
+    aPPs.push_back( *itPP );
+    if ( theHasAngles && itAngles != theAngles.end() )
+      aPPs.back().SetAngle( *itAngles );
   }
 
-  TNodeOfNodeListMap mapNewNodes;
+  TNodeOfNodeListMap   mapNewNodes;
   TElemOfVecOfNnlmiMap mapElemNewNodes;
-  TElemOfElemListMap newElemsMap;
+  TElemOfElemListMap   newElemsMap;
   TIDSortedElemSet::iterator itElem;
-  double aX, aY, aZ;
-  int aNb;
-  SMDSAbs_ElementType aTypeE;
   // source elements for each generated one
   SMESH_SequenceOfElemPtr srcElems, srcNodes;
 
   // 3. Center of rotation aV0
   gp_Pnt aV0 = theRefPoint;
-  gp_XYZ aGC;
-  if ( !theHasRefPoint ) {
-    aNb = 0;
-    aGC.SetCoord( 0.,0.,0. );
+  if ( !theHasRefPoint )
+  {
+    gp_XYZ aGC( 0.,0.,0. );
+    TIDSortedElemSet newNodes;
 
     itElem = theElements.begin();
     for ( ; itElem != theElements.end(); itElem++ ) {
@@ -5389,25 +5523,14 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet&  theElements,
 
       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
       while ( itN->more() ) {
-        const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( itN->next() );
-        aX = node->X();
-        aY = node->Y();
-        aZ = node->Z();
-
-        if ( mapNewNodes.find( node ) == mapNewNodes.end() ) {
-          list<const SMDS_MeshNode*> aLNx;
-          mapNewNodes[node] = aLNx;
-          //
-          gp_XYZ aXYZ( aX, aY, aZ );
-          aGC += aXYZ;
-          ++aNb;
-        }
+        const SMDS_MeshElement* node = itN->next();
+        if ( newNodes.insert( node ).second )
+          aGC += SMESH_TNodeXYZ( node );
       }
     }
-    aGC /= aNb;
+    aGC /= newNodes.size();
     aV0.SetXYZ( aGC );
   } // if (!theHasRefPoint) {
-  mapNewNodes.clear();
 
   // 4. Processing the elements
   SMESHDS_Mesh* aMesh = GetMeshDS();
@@ -5415,7 +5538,7 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet&  theElements,
   for ( itElem = theElements.begin(); itElem != theElements.end(); itElem++ ) {
     // check element type
     const SMDS_MeshElement* elem = *itElem;
-    aTypeE = elem->GetType();
+    SMDSAbs_ElementType   aTypeE = elem->GetType();
     if ( !elem || ( aTypeE != SMDSAbs_Face && aTypeE != SMDSAbs_Edge ) )
       continue;
 
@@ -5437,8 +5560,6 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet&  theElements,
         list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
 
         // make new nodes
-        aX = node->X();  aY = node->Y(); aZ = node->Z();
-
         Standard_Real aAngle1x, aAngleT1T0, aTolAng;
         gp_Pnt aP0x, aP1x, aPN0, aPN1, aV0x, aV1x;
         gp_Ax1 anAx1, anAxT1T0;
@@ -5447,17 +5568,17 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet&  theElements,
         aTolAng=1.e-4;
 
         aV0x = aV0;
-        aPN0.SetCoord(aX, aY, aZ);
+        aPN0 = SMESH_TNodeXYZ( node );
 
         const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0];
         aP0x = aPP0.Pnt();
         aDT0x= aPP0.Tangent();
         //cout<<"j = 0   PP: Pnt("<<aP0x.X()<<","<<aP0x.Y()<<","<<aP0x.Z()<<")"<<endl;
 
-        for ( j = 1; j < aNbTP; ++j ) {
+        for ( int j = 1; j < aNbTP; ++j ) {
           const SMESH_MeshEditor_PathPoint& aPP1 = aPPs[j];
-          aP1x = aPP1.Pnt();
-          aDT1x = aPP1.Tangent();
+          aP1x     = aPP1.Pnt();
+          aDT1x    = aPP1.Tangent();
           aAngle1x = aPP1.Angle();
 
           gp_Trsf aTrsf, aTrsfRot, aTrsfRotT1T0;
@@ -5501,10 +5622,7 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet&  theElements,
             srcNodes.Append( node );
             listNewNodes.push_back( newNode );
           }
-          aX = aPN1.X();
-          aY = aPN1.Y();
-          aZ = aPN1.Z();
-          const SMDS_MeshNode* newNode = aMesh->AddNode( aX, aY, aZ );
+          const SMDS_MeshNode* newNode = aMesh->AddNode( aPN1.X(), aPN1.Y(), aPN1.Z() );
           myLastCreatedNodes.Append(newNode);
           srcNodes.Append( node );
           listNewNodes.push_back( newNode );
@@ -8721,8 +8839,7 @@ void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements)
   }
 
   // replace given elements by linear ones
-  typedef SMDS_SetIterator<const SMDS_MeshElement*, TIDSortedElemSet::iterator> TSetIterator;
-  SMDS_ElemIteratorPtr elemIt( new TSetIterator( theElements.begin(), theElements.end() ));
+  SMDS_ElemIteratorPtr elemIt = elemSetIterator( theElements );
   removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 );
 
   // we need to convert remaining elements whose all medium nodes are in mediumNodeIDs
@@ -8774,8 +8891,7 @@ void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements)
       }
     }
   }
-  elemIt = SMDS_ElemIteratorPtr
-    (new TSetIterator( moreElemsToConvert.begin(), moreElemsToConvert.end() ));
+  elemIt = elemSetIterator( moreElemsToConvert );
   removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 );
 }
 
@@ -9460,6 +9576,95 @@ SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
   return SEW_OK;
 }
 
+//================================================================================
+/*!
+ * \brief Create elements equal (on same nodes) to given ones
+ *  \param [in] theElements - a set of elems to duplicate. If it is empty, all
+ *              elements of the uppest dimension are duplicated.
+ */
+//================================================================================
+
+void SMESH_MeshEditor::DoubleElements( const TIDSortedElemSet& theElements )
+{
+  CrearLastCreated();
+  SMESHDS_Mesh* mesh = GetMeshDS();
+
+  // get an element type and an iterator over elements
+
+  SMDSAbs_ElementType type;
+  SMDS_ElemIteratorPtr elemIt;
+  vector< const SMDS_MeshElement* > allElems;
+  if ( theElements.empty() )
+  {
+    if ( mesh->NbNodes() == 0 )
+      return;
+    // get most complex type
+    SMDSAbs_ElementType types[SMDSAbs_NbElementTypes] = {
+      SMDSAbs_Volume, SMDSAbs_Face, SMDSAbs_Edge,
+      SMDSAbs_0DElement, SMDSAbs_Ball, SMDSAbs_Node
+    };
+    for ( int i = 0; i < SMDSAbs_NbElementTypes; ++i )
+      if ( mesh->GetMeshInfo().NbElements( types[i] ))
+      {
+        type = types[i];
+        break;
+      }
+    // put all elements in the vector <allElems>
+    allElems.reserve( mesh->GetMeshInfo().NbElements( type ));
+    elemIt = mesh->elementsIterator( type );
+    while ( elemIt->more() )
+      allElems.push_back( elemIt->next());
+    elemIt = elemSetIterator( allElems );
+  }
+  else
+  {
+    type = (*theElements.begin())->GetType();
+    elemIt = elemSetIterator( theElements );
+  }
+
+  // duplicate elements
+
+  if ( type == SMDSAbs_Ball )
+  {
+    SMDS_UnstructuredGrid* vtkGrid = mesh->getGrid();
+    while ( elemIt->more() )
+    {
+      const SMDS_MeshElement* elem = elemIt->next();
+      if ( elem->GetType() != SMDSAbs_Ball )
+        continue;
+      if (( elem = mesh->AddBall( elem->GetNode(0),
+                                  vtkGrid->GetBallDiameter( elem->getVtkId() ))))
+        myLastCreatedElems.Append( elem );
+    }
+  }
+  else
+  {
+    vector< const SMDS_MeshNode* > nodes;
+    while ( elemIt->more() )
+    {
+      const SMDS_MeshElement* elem = elemIt->next();
+      if ( elem->GetType() != type )
+        continue;
+
+      nodes.assign( elem->begin_nodes(), elem->end_nodes() );
+
+      if ( type == SMDSAbs_Volume  && elem->GetVtkType() == VTK_POLYHEDRON )
+      {
+        std::vector<int> quantities =
+          static_cast< const SMDS_VtkVolume* >( elem )->GetQuantities();
+        elem = mesh->AddPolyhedralVolume( nodes, quantities );
+      }
+      else
+      {
+        AddElement( nodes, type, elem->IsPoly() );
+        elem = 0; // myLastCreatedElems is already filled
+      }
+      if ( elem )
+        myLastCreatedElems.Append( elem );
+    }
+  }
+}
+
 //================================================================================
 /*!
   \brief Creates a hole in a mesh by doubling the nodes of some particular elements
@@ -10764,7 +10969,8 @@ void SMESH_MeshEditor::CreateHoleSkin(double radius,
           double z = nodesCoords[i++];
           gp_Pnt p = gp_Pnt(x, y ,z);
           gpnts.push_back(p);
-          MESSAGE("TopoDS_Vertex " << k++ << " " << p.X() << " " << p.Y() << " " << p.Z());
+          MESSAGE("TopoDS_Vertex " << k << " " << p.X() << " " << p.Y() << " " << p.Z());
+          k++;
         }
     }
   else // --- no group, no coordinates : use the vertices of the geom shape provided, and radius
@@ -11262,10 +11468,8 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
   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() ));
+  if (elements.empty()) eIt = aMesh->elementsIterator(elemType);
+  else                  eIt = elemSetIterator( elements );
 
   while (eIt->more())
   {
@@ -11459,10 +11663,8 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
   // -----------------------
   if ( toCopyElements && targetMesh != myMesh )
   {
-    if (elements.empty())
-      eIt = aMesh->elementsIterator(elemType);
-    else
-      eIt = SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() ));
+    if (elements.empty()) eIt = aMesh->elementsIterator(elemType);
+    else                  eIt = elemSetIterator( elements );
     while (eIt->more())
     {
       const SMDS_MeshElement* elem = eIt->next();