Salome HOME
Copyright update 2022
[modules/smesh.git] / src / SMESH / SMESH_MeshEditor.cxx
index 22695bac0567cc812b241097973dc6ae7a5bc028..be874359b73c68b2ab7891c3acb3c09104a2984c 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2021  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2022  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
@@ -412,8 +412,8 @@ SMDS_MeshElement* SMESH_MeshEditor::AddElement(const vector<smIdType> & nodeIDs,
 //           Modify a compute state of sub-meshes which become empty
 //=======================================================================
 
-smIdType SMESH_MeshEditor::Remove (const list< smIdType >& theIDs,
-                              const bool         isNodes )
+smIdType SMESH_MeshEditor::Remove (const std::list< smIdType >& theIDs,
+                                   const bool                   isNodes )
 {
   ClearLastCreated();
 
@@ -474,10 +474,87 @@ smIdType SMESH_MeshEditor::Remove (const list< smIdType >& theIDs,
   return removed;
 }
 
+//================================================================================
+/*!
+ * \brief Remove a node and fill a hole appeared, by changing surrounding faces
+ */
+//================================================================================
+
+void SMESH_MeshEditor::RemoveNodeWithReconnection( const SMDS_MeshNode* node )
+{
+  if ( ! node )
+    return;
+
+  if ( node->NbInverseElements( SMDSAbs_Volume ) > 0 )
+    throw SALOME_Exception( "RemoveNodeWithReconnection() applies to 2D mesh only" );
+
+  // check that only triangles surround the node
+  for ( SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator( SMDSAbs_Face ); fIt->more(); )
+  {
+    const SMDS_MeshElement* face = fIt->next();
+    if ( face->GetGeomType() != SMDSGeom_TRIANGLE )
+      throw SALOME_Exception( "RemoveNodeWithReconnection() applies to triangle mesh only" );
+    if ( face->IsQuadratic() )
+      throw SALOME_Exception( "RemoveNodeWithReconnection() applies to linear mesh only" );
+  }
+
+  std::vector< const SMDS_MeshNode*> neighbours(2);
+  SMESH_MeshAlgos::IsOn2DBoundary( node, & neighbours );
+
+  bool toRemove = ( neighbours.size() > 2 ); // non-manifold ==> just remove
+
+  // if ( neighbours.size() == 2 ) // on boundary
+  // {
+  //   // check if theNode and neighbours are on a line
+  //   gp_Pnt pN = SMESH_NodeXYZ( node );
+  //   gp_Pnt p0 = SMESH_NodeXYZ( neighbours[0] );
+  //   gp_Pnt p1 = SMESH_NodeXYZ( neighbours[1] );
+  //   double dist01 = p0.Distance( p1 );
+  //   double    tol = 0.01 * dist01;
+  //   double  distN = ( gp_Vec( p0, p1 ) ^ gp_Vec( p0, pN )).Magnitude() / dist01;
+  //   bool   onLine = distN < tol;
+  //   toRemove = !onLine;
+  // }
+
+  if ( neighbours.empty() ) // not on boundary
+  {
+    TIDSortedElemSet linkedNodes;
+    GetLinkedNodes( node, linkedNodes, SMDSAbs_Face );
+    for ( const SMDS_MeshElement* e : linkedNodes ) neighbours.push_back( cast2Node( e ));
+    if ( neighbours.empty() )
+      toRemove = true;
+  }
+
+  if ( toRemove )
+  {
+    this->Remove( std::list< smIdType >( 1, node->GetID() ), /*isNode=*/true );
+    return;
+  }
+
+  // choose a node to replace by
+  const SMDS_MeshNode* nToReplace = nullptr;
+  SMESH_NodeXYZ           nodeXYZ = node;
+  double                  minDist = Precision::Infinite();
+  for ( const SMDS_MeshNode* n : neighbours )
+  {
+    double dist = nodeXYZ.SquareDistance( n );
+    if ( dist < minDist )
+    {
+      minDist = dist;
+      nToReplace = n;
+    }
+  }
+
+  // remove node + replace by nToReplace
+  std::list< const SMDS_MeshNode* > nodeGroup = { nToReplace, node };
+  TListOfListOfNodes nodesToMerge( 1, nodeGroup );
+  this->MergeNodes( nodesToMerge );
+}
+
 //================================================================================
 /*!
  * \brief Create 0D elements on all nodes of the given object.
- *  \param elements - Elements on whose nodes to create 0D elements; if empty, 
+ *  \param elements - Elements on whose nodes to create 0D elements; if empty,
  *                    the all mesh is treated
  *  \param all0DElems - returns all 0D elements found or created on nodes of \a elements
  *  \param duplicateElements - to add one more 0D element to a node or not
@@ -947,24 +1024,24 @@ bool getQuadrangleNodes(const SMDS_MeshNode *    theQuadNodes [],
 {
   if( tr1->NbNodes() != tr2->NbNodes() )
     return false;
+
   // find the 4-th node to insert into tr1
   const SMDS_MeshNode* n4 = 0;
   SMDS_ElemIteratorPtr it = tr2->nodesIterator();
-  int i=0;
-  while ( !n4 && i<3 ) {
+  for ( int i = 0; !n4 && i < 3; ++i )
+  {
     const SMDS_MeshNode * n = cast2Node( it->next() );
-    i++;
     bool isDiag = ( n == theNode1 || n == theNode2 );
     if ( !isDiag )
       n4 = n;
   }
+
   // Make an array of nodes to be in a quadrangle
   int iNode = 0, iFirstDiag = -1;
   it = tr1->nodesIterator();
-  i=0;
-  while ( i<3 ) {
+  for ( int i = 0; i < 3; ++i )
+  {
     const SMDS_MeshNode * n = cast2Node( it->next() );
-    i++;
     bool isDiag = ( n == theNode1 || n == theNode2 );
     if ( isDiag ) {
       if ( iFirstDiag < 0 )
@@ -1079,6 +1156,210 @@ bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
   return true;
 }
 
+//=======================================================================
+//function : SplitEdge
+//purpose  : Replace each triangle bound by theNode1-theNode2 segment with
+//           two triangles by connecting a node made on the link with a node opposite to the link.
+//=======================================================================
+
+void SMESH_MeshEditor::SplitEdge (const SMDS_MeshNode * theNode1,
+                                  const SMDS_MeshNode * theNode2,
+                                  double                thePosition)
+{
+  ClearLastCreated();
+
+  SMESHDS_Mesh * mesh = GetMeshDS();
+
+  // Get triangles and segments to divide
+
+  std::vector<const SMDS_MeshNode *> diagNodes = { theNode1, theNode2 };
+  std::vector<const SMDS_MeshElement *> foundElems;
+  if ( !mesh->GetElementsByNodes( diagNodes, foundElems ) || foundElems.empty() )
+    throw SALOME_Exception( SMESH_Comment("No triangle is bound by the edge ")
+                            << theNode1->GetID() << " - " << theNode2->GetID());
+
+  SMESH_MesherHelper helper( *GetMesh() );
+
+  for ( const SMDS_MeshElement * elem : foundElems )
+  {
+    SMDSAbs_ElementType type = elem->GetType();
+    switch ( type ) {
+    case SMDSAbs_Volume:
+      throw SALOME_Exception( "Can't split an edge of a volume");
+      break;
+
+    case SMDSAbs_Face:
+      if ( elem->GetGeomType() != SMDSGeom_TRIANGLE )
+        throw SALOME_Exception( "Can't split an edge of a face of type other than triangle");
+      if ( elem->IsQuadratic() )
+      {
+        helper.SetIsQuadratic( true );
+        helper.AddTLinks( static_cast< const SMDS_MeshFace*>( elem ));
+        helper.SetIsBiQuadratic( elem->GetEntityType() == SMDSEntity_BiQuad_Triangle );
+      }
+      break;
+
+    case SMDSAbs_Edge:
+      if ( elem->IsQuadratic() )
+      {
+        helper.SetIsQuadratic( true );
+        helper.AddTLinks( static_cast< const SMDS_MeshEdge*>( elem ));
+      }
+      break;
+    default:;
+    }
+  }
+
+  // Make a new node
+
+  const SMDS_MeshNode* nodeOnLink = helper.GetMediumNode( theNode1, theNode2,/*force3d=*/false );
+
+  gp_Pnt newNodeXYZ = ( SMESH_NodeXYZ( theNode1 ) * ( 1 - thePosition ) +
+                        SMESH_NodeXYZ( theNode2 ) * thePosition );
+
+  const TopoDS_Shape& S = mesh->IndexToShape( nodeOnLink->GetShapeID() );
+  if ( !S.IsNull() && S.ShapeType() == TopAbs_FACE ) // find newNodeXYZ by UV on FACE
+  {
+    Handle(ShapeAnalysis_Surface) surface = helper.GetSurface( TopoDS::Face( S ));
+    double  tol = 100 * helper.MaxTolerance( S );
+    gp_Pnt2d uv = surface->ValueOfUV( newNodeXYZ, tol );
+    if ( surface->Gap() < SMESH_NodeXYZ( theNode1 ).Distance( theNode2 ))
+    {
+      newNodeXYZ = surface->Value( uv );
+      if ( SMDS_FacePositionPtr nPos = nodeOnLink->GetPosition())
+        nPos->SetParameters( uv.X(), uv.Y() );
+    }
+  }
+  if ( !S.IsNull() && S.ShapeType() == TopAbs_EDGE ) // find newNodeXYZ by param on EDGE
+  {
+    mesh->MoveNode( nodeOnLink, newNodeXYZ.X(), newNodeXYZ.Y(), newNodeXYZ.Z() );
+    double u = Precision::Infinite(), tol = 100 * helper.MaxTolerance( S ), distXYZ[4];
+    helper.ToFixNodeParameters( true );
+    if ( helper.CheckNodeU( TopoDS::Edge( S ), nodeOnLink, u, tol, /*force3D=*/false, distXYZ ))
+      newNodeXYZ.SetCoord( distXYZ[1], distXYZ[2], distXYZ[3] );
+  }
+  mesh->MoveNode( nodeOnLink, newNodeXYZ.X(), newNodeXYZ.Y(), newNodeXYZ.Z() );
+
+  // Split triangles and segments
+
+  std::vector<const SMDS_MeshNode *> nodes( 7 );
+  for ( const SMDS_MeshElement * elem : foundElems )
+  {
+    nodes.assign( elem->begin_nodes(), elem->end_nodes() );
+    nodes.resize( elem->NbCornerNodes() + 1 );
+    nodes.back() = nodes[0];
+
+    smIdType id = elem->GetID();
+    int shapeID = elem->GetShapeID();
+
+    const SMDS_MeshNode* centralNode = nullptr;
+    if ( elem->GetEntityType() == SMDSEntity_BiQuad_Triangle )
+      centralNode = elem->GetNode( 6 );
+
+    mesh->RemoveFreeElement( elem, /*sm=*/0, /*fromGroups=*/false );
+    if ( centralNode )
+      mesh->RemoveFreeNode( centralNode, /*sm=*/0, /*fromGroups=*/true );
+
+    for ( size_t i = 1; i < nodes.size(); ++i )
+    {
+      const SMDS_MeshNode* n1 = nodes[i-1];
+      const SMDS_MeshNode* n2 = nodes[i];
+      const SMDS_MeshElement* newElem;
+      if ( nodes.size() == 4 ) //    triangle
+      {
+        bool isDiag1 = ( n1 == theNode1 || n1 == theNode2 );
+        bool isDiag2 = ( n2 == theNode1 || n2 == theNode2 );
+        if ( isDiag1 && isDiag2 )
+          continue;
+
+        newElem = helper.AddFace( n1, n2, nodeOnLink, id );
+      }
+      else //    segment
+      {
+        newElem = helper.AddEdge( n1, nodeOnLink, id );
+      }
+      myLastCreatedElems.push_back( newElem );
+      AddToSameGroups( newElem, elem, mesh );
+      if ( shapeID )
+        mesh->SetMeshElementOnShape( newElem, shapeID );
+      id = 0;
+    }
+  }
+  return;
+}
+
+//=======================================================================
+//function : SplitFace
+//purpose  : Split a face into triangles each formed by two nodes of the 
+//           face and a new node added at the given coordinates.
+//=======================================================================
+
+void SMESH_MeshEditor::SplitFace (const SMDS_MeshElement * theFace,
+                                  double                   theX,
+                                  double                   theY,
+                                  double                   theZ )
+{
+  ClearLastCreated();
+
+  if ( !theFace )
+    throw SALOME_Exception("Null face given");
+  if ( theFace->GetType() != SMDSAbs_Face )
+    throw SALOME_Exception("Not a face given");
+
+  SMESHDS_Mesh * mesh = GetMeshDS();
+
+  SMESH_MesherHelper helper( *GetMesh() );
+  if ( theFace->IsQuadratic() )
+  {
+    helper.SetIsQuadratic( true );
+    helper.AddTLinks( static_cast< const SMDS_MeshFace*>( theFace ));
+  }
+  const TopoDS_Shape& shape = mesh->IndexToShape( theFace->GetShapeID() );
+  helper.SetSubShape( shape );
+  helper.SetElementsOnShape( true );
+
+  // Make a new node
+
+  const SMDS_MeshNode* centralNode = nullptr;
+  if (      theFace->GetEntityType() == SMDSEntity_BiQuad_Triangle )
+    centralNode = theFace->GetNode( 6 );
+  else if ( theFace->GetEntityType() == SMDSEntity_BiQuad_Quadrangle )
+    centralNode = theFace->GetNode( 8 );
+
+  if ( centralNode )
+  {
+    helper.SetIsBiQuadratic( true );
+    mesh->MoveNode( centralNode, theX, theY, theZ );
+  }
+  else
+    centralNode = helper.AddNode( theX, theY, theZ );
+
+
+  // Split theFace
+
+  std::vector<const SMDS_MeshNode *> nodes( theFace->NbNodes() + 1 );
+  nodes.assign( theFace->begin_nodes(), theFace->end_nodes() );
+  nodes.resize( theFace->NbCornerNodes() + 1 );
+  nodes.back() = nodes[0];
+
+  smIdType id = theFace->GetID();
+  int shapeID = theFace->GetShapeID();
+
+  mesh->RemoveFreeElement( theFace, /*sm=*/0, /*fromGroups=*/false );
+
+  for ( size_t i = 1; i < nodes.size(); ++i )
+  {
+    const SMDS_MeshElement* newElem = helper.AddFace( nodes[i-1], nodes[i], centralNode, id );
+
+    myLastCreatedElems.push_back( newElem );
+    AddToSameGroups( newElem, theFace, mesh );
+    if ( shapeID )
+      mesh->SetMeshElementOnShape( newElem, shapeID );
+    id = 0;
+  }
+  return;
+}
+
 //=======================================================================
 //function : Reorient
 //purpose  : Reverse theElement orientation
@@ -1161,69 +1442,88 @@ bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
 /*!
  * \brief Reorient faces.
  * \param theFaces - the faces to reorient. If empty the whole mesh is meant
- * \param theDirection - desired direction of normal of \a theFace
- * \param theFace - one of \a theFaces that should be oriented according to
- *        \a theDirection and whose orientation defines orientation of other faces
+ * \param theDirection - desired direction of normal of \a theRefFaces.
+ *        It can be (0,0,0) in order to keep orientation of \a theRefFaces.
+ * \param theRefFaces - correctly oriented faces whose orientation defines
+ *        orientation of other faces.
  * \return number of reoriented faces.
  */
 //================================================================================
 
-int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
-                                  const gp_Dir&            theDirection,
-                                  const SMDS_MeshElement * theFace)
+int SMESH_MeshEditor::Reorient2D( TIDSortedElemSet &  theFaces,
+                                  const gp_Vec&       theDirection,
+                                  TIDSortedElemSet &  theRefFaces,
+                                  bool                theAllowNonManifold )
 {
   int nbReori = 0;
-  if ( !theFace || theFace->GetType() != SMDSAbs_Face ) return nbReori;
 
   if ( theFaces.empty() )
   {
-    SMDS_FaceIteratorPtr fIt = GetMeshDS()->facesIterator(/*idInceasingOrder=true*/);
+    SMDS_FaceIteratorPtr fIt = GetMeshDS()->facesIterator();
     while ( fIt->more() )
       theFaces.insert( theFaces.end(), fIt->next() );
+
+    if ( theFaces.empty() )
+      return nbReori;
   }
 
-  // orient theFace according to theDirection
-  gp_XYZ normal;
-  SMESH_MeshAlgos::FaceNormal( theFace, normal, /*normalized=*/false );
-  if ( normal * theDirection.XYZ() < 0 )
-    nbReori += Reorient( theFace );
+  // orient theRefFaces according to theDirection
+  if ( theDirection.X() != 0 || theDirection.Y() != 0 || theDirection.Z() != 0 )
+    for ( const SMDS_MeshElement* refFace : theRefFaces )
+    {
+      gp_XYZ normal;
+      SMESH_MeshAlgos::FaceNormal( refFace, normal, /*normalized=*/false );
+      if ( normal * theDirection.XYZ() < 0 )
+        nbReori += Reorient( refFace );
+    }
 
-  // Orient other faces
+  // mark reference faces
+  GetMeshDS()->SetAllCellsNotMarked();
+  for ( const SMDS_MeshElement* refFace : theRefFaces )
+    refFace->setIsMarked( true );
 
-  set< const SMDS_MeshElement* > startFaces, visitedFaces;
-  TIDSortedElemSet avoidSet;
-  set< SMESH_TLink > checkedLinks;
-  pair< set< SMESH_TLink >::iterator, bool > linkIt_isNew;
+  // erase reference faces from theFaces
+  for ( TIDSortedElemSet::iterator fIt = theFaces.begin(); fIt != theFaces.end(); )
+    if ( (*fIt)->isMarked() )
+      fIt = theFaces.erase( fIt );
+    else
+      ++fIt;
 
-  if ( theFaces.size() > 1 )// leave 1 face to prevent finding not selected faces
-    theFaces.erase( theFace );
-  startFaces.insert( theFace );
+  if ( theRefFaces.empty() )
+  {
+    theRefFaces.insert( *theFaces.begin() );
+    theFaces.erase( theFaces.begin() );
+  }
+
+  // Orient theFaces
+
+  // if ( theFaces.size() > 1 )// leave 1 face to prevent finding not selected faces
+  //   theFaces.erase( theFace );
 
   int nodeInd1, nodeInd2;
-  const SMDS_MeshElement*           otherFace;
+  const SMDS_MeshElement*           refFace, *otherFace;
   vector< const SMDS_MeshElement* > facesNearLink;
   vector< std::pair< int, int > >   nodeIndsOfFace;
+  TIDSortedElemSet                  avoidSet, emptySet;
+  NCollection_Map< SMESH_TLink, SMESH_TLink > checkedLinks;
 
-  set< const SMDS_MeshElement* >::iterator startFace = startFaces.begin();
-  while ( !startFaces.empty() )
+  while ( !theRefFaces.empty() )
   {
-    startFace = startFaces.begin();
-    theFace = *startFace;
-    startFaces.erase( startFace );
-    if ( !visitedFaces.insert( theFace ).second )
-      continue;
+    auto refFaceIt = theRefFaces.begin();
+    refFace = *refFaceIt;
+    theRefFaces.erase( refFaceIt );
 
     avoidSet.clear();
-    avoidSet.insert(theFace);
+    avoidSet.insert( refFace );
 
-    NLink link( theFace->GetNode( 0 ), (SMDS_MeshNode *) 0 );
+    NLink link( refFace->GetNode( 0 ), nullptr );
 
-    const int nbNodes = theFace->NbCornerNodes();
-    for ( int i = 0; i < nbNodes; ++i ) // loop on links of theFace
+    const int nbNodes = refFace->NbCornerNodes();
+    for ( int i = 0; i < nbNodes; ++i ) // loop on links of refFace
     {
-      link.second = theFace->GetNode(( i+1 ) % nbNodes );
-      linkIt_isNew = checkedLinks.insert( link );
-      if ( !linkIt_isNew.second )
+      link.second = refFace->GetNode(( i+1 ) % nbNodes );
+      bool isLinkVisited = checkedLinks.Contains( link );
+      if ( isLinkVisited )
       {
         // link has already been checked and won't be encountered more
         // if the group (theFaces) is manifold
@@ -1231,28 +1531,41 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
       }
       else
       {
+        checkedLinks.Add( link );
+
         facesNearLink.clear();
         nodeIndsOfFace.clear();
+        TIDSortedElemSet::iterator objFaceIt = theFaces.end();
+
         while (( otherFace = SMESH_MeshAlgos::FindFaceInSet( link.first, link.second,
-                                                             theFaces, avoidSet,
+                                                             emptySet, avoidSet,
                                                              &nodeInd1, &nodeInd2 )))
-          if ( otherFace != theFace)
+        {
+          if (( otherFace->isMarked() ) || // ref face
+              (( objFaceIt = theFaces.find( otherFace )) != theFaces.end() )) // object face
           {
             facesNearLink.push_back( otherFace );
             nodeIndsOfFace.push_back( make_pair( nodeInd1, nodeInd2 ));
-            avoidSet.insert( otherFace );
           }
+          avoidSet.insert( otherFace );
+        }
         if ( facesNearLink.size() > 1 )
         {
           // NON-MANIFOLD mesh shell !
-          // select a face most co-directed with theFace,
+          if ( !theAllowNonManifold )
+          {
+            throw SALOME_Exception("Non-manifold topology of groups");
+          }
+          // select a face most co-directed with refFace,
           // other faces won't be visited this time
           gp_XYZ NF, NOF;
-          SMESH_MeshAlgos::FaceNormal( theFace, NF, /*normalized=*/false );
+          SMESH_MeshAlgos::FaceNormal( refFace, NF, /*normalized=*/false );
           double proj, maxProj = -1;
-          for ( size_t i = 0; i < facesNearLink.size(); ++i ) {
+          for ( size_t i = 0; i < facesNearLink.size(); ++i )
+          {
             SMESH_MeshAlgos::FaceNormal( facesNearLink[i], NOF, /*normalized=*/false );
-            if (( proj = Abs( NF * NOF )) > maxProj ) {
+            if (( proj = Abs( NF * NOF )) > maxProj )
+            {
               maxProj = proj;
               otherFace = facesNearLink[i];
               nodeInd1  = nodeIndsOfFace[i].first;
@@ -1260,9 +1573,9 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
             }
           }
           // not to visit rejected faces
-          for ( size_t i = 0; i < facesNearLink.size(); ++i )
-            if ( facesNearLink[i] != otherFace && theFaces.size() > 1 )
-              visitedFaces.insert( facesNearLink[i] );
+          // for ( size_t i = 0; i < facesNearLink.size(); ++i )
+          //   if ( facesNearLink[i] != otherFace && theFaces.size() > 1 )
+          //     visitedFaces.insert( facesNearLink[i] );
         }
         else if ( facesNearLink.size() == 1 )
         {
@@ -1270,20 +1583,36 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
           nodeInd1  = nodeIndsOfFace.back().first;
           nodeInd2  = nodeIndsOfFace.back().second;
         }
-        if ( otherFace && otherFace != theFace)
+        if ( otherFace )
         {
-          // link must be reverse in otherFace if orientation to otherFace
-          // is same as that of theFace
-          if ( abs(nodeInd2-nodeInd1) == 1 ? nodeInd2 > nodeInd1 : nodeInd1 > nodeInd2 )
+          // link must be reverse in otherFace if orientation of otherFace
+          // is same as that of refFace
+          if ( abs( nodeInd2 - nodeInd1 ) == 1 ? nodeInd2 > nodeInd1 : nodeInd1 > nodeInd2 )
           {
+            if ( otherFace->isMarked() )
+              throw SALOME_Exception("Different orientation of reference faces");
             nbReori += Reorient( otherFace );
           }
-          startFaces.insert( otherFace );
+          if ( !otherFace->isMarked() )
+          {
+            theRefFaces.insert( otherFace );
+            if ( objFaceIt != theFaces.end() )
+              theFaces.erase( objFaceIt );
+          }
         }
       }
-      std::swap( link.first, link.second ); // reverse the link
+      link.first = link.second; // reverse the link
+
+    } // loop on links of refFace
+
+    if ( theRefFaces.empty() && !theFaces.empty() )
+    {
+      theRefFaces.insert( *theFaces.begin() );
+      theFaces.erase( theFaces.begin() );
     }
-  }
+
+  } // while ( !theRefFaces.empty() )
+
   return nbReori;
 }
 
@@ -6959,9 +7288,19 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes,
 
     for ( size_t i = 0; i < newElemDefs.size(); ++i )
     {
-      if ( i > 0 || !mesh->ChangeElementNodes( elem,
-                                               & newElemDefs[i].myNodes[0],
-                                               newElemDefs[i].myNodes.size() ))
+      bool elemChanged = false;
+      if ( i == 0 )
+      {
+        if ( elem->GetGeomType() == SMDSGeom_POLYHEDRA )
+          elemChanged = mesh->ChangePolyhedronNodes( elem,
+                                                     newElemDefs[i].myNodes,
+                                                     newElemDefs[i].myPolyhedQuantities );
+        else
+          elemChanged = mesh->ChangeElementNodes( elem,
+                                                  & newElemDefs[i].myNodes[0],
+                                                  newElemDefs[i].myNodes.size() );
+      }
+      if ( i > 0 || !elemChanged )
       {
         if ( i == 0 )
         {
@@ -7128,6 +7467,7 @@ bool SMESH_MeshEditor::applyMerge( const SMDS_MeshElement* elem,
         // each face has to be analyzed in order to check volume validity
         if ( const SMDS_MeshVolume* aPolyedre = SMDS_Mesh::DownCast< SMDS_MeshVolume >( elem ))
         {
+          toRemove = false;
           int nbFaces = aPolyedre->NbFaces();
 
           vector<const SMDS_MeshNode *>& poly_nodes = newElemDefs[0].myNodes;
@@ -10781,7 +11121,23 @@ bool SMESH_MeshEditor::doubleNodes(SMESHDS_Mesh*           theMeshDS,
     if ( theIsDoubleElem )
       AddElement( newNodes, elemType.Init( anElem, /*basicOnly=*/false ));
     else
-      theMeshDS->ChangeElementNodes( anElem, &newNodes[ 0 ], newNodes.size() );
+    {
+      // change element nodes
+      const SMDSAbs_EntityType geomType = anElem->GetEntityType();
+      if ( geomType == SMDSEntity_Polyhedra )
+      {
+        // special treatment for polyhedron
+        const SMDS_MeshVolume* aPolyhedron = SMDS_Mesh::DownCast< SMDS_MeshVolume >( anElem );
+        if (!aPolyhedron) {
+          MESSAGE("Warning: bad volumic element");
+          return false;
+        }
+        theMeshDS->ChangePolyhedronNodes( anElem, newNodes, aPolyhedron->GetQuantities() );
+      }
+      else
+        // standard entity type
+        theMeshDS->ChangeElementNodes( anElem, &newNodes[ 0 ], newNodes.size() );
+    }
 
     res = true;
   }
@@ -11107,16 +11463,9 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   //     build the list of nodes shared by 2 or more domains, with their domain indexes
 
   std::map<DownIdType, std::map<int,int>, DownIdCompare> faceDomains; // face --> (id domain --> id volume)
-  std::map<int,int>celldom; // cell vtkId --> domain
+  std::map<int,int> celldom; // cell vtkId --> domain
   std::map<DownIdType, std::map<int,int>, DownIdCompare> cellDomains;  // oldNode --> (id domain --> id cell)
   std::map<int, std::map<int,int> > nodeDomains; // oldId -->  (domainId --> newId)
-  faceDomains.clear();
-  celldom.clear();
-  cellDomains.clear();
-  nodeDomains.clear();
-  std::map<int,int> emptyMap;
-  std::set<int> emptySet;
-  emptyMap.clear();
 
   //MESSAGE(".. Number of domains :"<<theElems.size());
 
@@ -11226,7 +11575,6 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       DownIdType face = itface->first;
       //MESSAGE(" --- face " << face.cellId);
       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)
@@ -11281,7 +11629,6 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       DownIdType face = itface->first;
       //MESSAGE(" --- face " << face.cellId);
       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)
@@ -11342,7 +11689,6 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       DownIdType face = itface->first;
       //MESSAGE(" --- face " << face.cellId);
       std::set<int> oldNodes;
-      oldNodes.clear();
       grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
       int nbMultipleNodes = 0;
       std::set<int>::iterator itn = oldNodes.begin();
@@ -11467,7 +11813,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
 
   //MESSAGE(".. Creation of elements: simple junction");
-  if (createJointElems)
+  if ( createJointElems )
   {
     string joints2DName = "joints2D";
     mapOfJunctionGroups[joints2DName] = this->myMesh->AddGroup(SMDSAbs_Face, joints2DName.c_str());
@@ -11482,15 +11828,14 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       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> domvol = itface->second;
+      std::map<int, int>          domvol = itface->second;
       std::map<int, int>::iterator itdom = domvol.begin();
-      int dom1 = itdom->first;
+      int     dom1 = itdom->first;
       int vtkVolId = itdom->second;
       itdom++;
-      int dom2 = itdom->first;
+      int           dom2 = itdom->first;
       SMDS_MeshCell *vol = grid->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains,
                                                        nodeQuadDomains);
       stringstream grpname;
@@ -11547,7 +11892,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
     std::map<std::vector<int>, std::vector<int> >::iterator ite = edgesMultiDomains.begin();
     for (; ite != edgesMultiDomains.end(); ++ite)
     {
-      vector<int> nodes = ite->first;
+      vector<int>    nodes = ite->first;
       vector<int> orderDom = ite->second;
       vector<vtkIdType> orderedNodes;
       if (nodes.size() == 2)
@@ -11584,10 +11929,9 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
 
   std::map<DownIdType, std::map<int,int>, DownIdCompare> faceOrEdgeDom; // cellToModify --> (id domain --> id cell)
   std::map<int,int> feDom; // vtk id of cell to modify --> id domain
-  faceOrEdgeDom.clear();
-  feDom.clear();
 
   //MESSAGE(".. Modification of elements");
+  SMDSAbs_ElementType domainType = (*theElems[0].begin())->GetType();
   for (int idomain = idom0; idomain < nbDomains; idomain++)
   {
     std::map<int, std::map<int, int> >::const_iterator itnod = nodeDomains.begin();
@@ -11605,13 +11949,29 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
           continue; // new cells: not to be modified
         DownIdType aCell(downId, vtkType);
         int volParents[1000];
-        int nbvol = grid->GetParentVolumes(volParents, vtkId);
+        int nbvol = 0;
+        nbvol = grid->GetParentVolumes(volParents, vtkId);
+        if ( domainType == SMDSAbs_Volume )
+        {
+          nbvol = grid->GetParentVolumes(volParents, vtkId);
+        }
+        else // domainType == SMDSAbs_Face
+        {
+          const int            nbFaces = grid->getDownArray(vtkType)->getNumberOfUpCells(downId);
+          const int           *upCells = grid->getDownArray(vtkType)->getUpCells(downId);
+          const unsigned char* upTypes = grid->getDownArray(vtkType)->getUpTypes(downId);
+          for (int i=0; i< nbFaces; i++)
+          {
+            int vtkFaceId = grid->getDownArray( upTypes[i] )->getVtkCellId(upCells[i]);
+            if (vtkFaceId >= 0)
+              volParents[nbvol++] = vtkFaceId;
+          }
+        }
         for (int j = 0; j < nbvol; j++)
           if (celldom.count(volParents[j]) && (celldom[volParents[j]] == idomain))
             if (!feDom.count(vtkId))
             {
               feDom[vtkId] = idomain;
-              faceOrEdgeDom[aCell] = emptyMap;
               faceOrEdgeDom[aCell][idomain] = vtkId; // affect face or edge to the first domain only
               //MESSAGE("affect cell " << this->GetMeshDS()->FromVtkToSmds(vtkId) << " domain " << idomain
               //        << " type " << vtkType << " downId " << downId);
@@ -11634,7 +11994,6 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       DownIdType face = itface->first;
       std::set<int> oldNodes;
       std::set<int>::iterator itn;
-      oldNodes.clear();
       grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
       //MESSAGE("examine cell, downId " << face.cellId << " type " << int(face.cellType));
       std::map<int, int> localClonedNodeIds;
@@ -11701,10 +12060,7 @@ bool SMESH_MeshEditor::CreateFlatElementsOnFacesGroups(const std::vector<TIDSort
 
   std::map<const SMDS_MeshNode*, const SMDS_MeshNode*> clonedNodes;
   std::map<const SMDS_MeshNode*, const SMDS_MeshNode*> intermediateNodes;
-  clonedNodes.clear();
-  intermediateNodes.clear();
   std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
-  mapOfJunctionGroups.clear();
 
   for ( size_t idom = 0; idom < theElems.size(); idom++ )
   {
@@ -11951,7 +12307,6 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
   std::set<int> setOfVolToCheck;
 
   std::vector<gp_Pnt> gpnts;
-  gpnts.clear();
 
   if (isNodeGroup) // --- a group of nodes is provided : find all the volumes using one or more of this nodes
   {
@@ -12038,7 +12393,6 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
     //     Fill the group of inside volumes
 
     std::map<int, double> mapOfNodeDistance2;
-    mapOfNodeDistance2.clear();
     std::set<int> setOfOutsideVol;
     while (!setOfVolToCheck.empty())
     {
@@ -12215,9 +12569,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
   //     project polylines on subshapes, and partition, to get geom faces
 
   std::map<int, std::set<int> > shapeIdToVtkIdSet; // shapeId --> set of vtkId on skin
-  std::set<int> emptySet;
-  emptySet.clear();
-  std::set<int> shapeIds;
+  std::set<int>                 shapeIds;
 
   SMDS_ElemIteratorPtr itelem = sgrps->GetElements();
   while (itelem->more())
@@ -12227,7 +12579,6 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
     int   vtkId = elem->GetVtkID();
     if (!shapeIdToVtkIdSet.count(shapeId))
     {
-      shapeIdToVtkIdSet[shapeId] = emptySet;
       shapeIds.insert(shapeId);
     }
     shapeIdToVtkIdSet[shapeId].insert(vtkId);
@@ -12235,7 +12586,6 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
 
   std::map<int, std::set<DownIdType, DownIdCompare> > shapeIdToEdges; // shapeId --> set of downward edges
   std::set<DownIdType, DownIdCompare> emptyEdges;
-  emptyEdges.clear();
 
   std::map<int, std::set<int> >::iterator itShape =  shapeIdToVtkIdSet.begin();
   for (; itShape != shapeIdToVtkIdSet.end(); ++itShape)
@@ -12278,7 +12628,6 @@ void SMESH_MeshEditor::CreateHoleSkin(double                          radius,
     }
 
     std::list<int> order;
-    order.clear();
     if (nodesEdges.size() > 0)
     {
       order.push_back(nodesEdges[0]); //MESSAGE("       --- back " << order.back()+1); // SMDS id = VTK id + 1;