Salome HOME
Regression of XSMESH_TEST/SMESHCOMMON/SMESH_TEST/Grids/smesh/bugs12/M6
[modules/smesh.git] / src / SMESH / SMESH_MeshEditor.cxx
index 9a821ab80c4e221cc4b32ccaad50fca18265b09d..f5bc661fbe1db0e0597e36aff91ac17280bab22c 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2013  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
@@ -97,6 +97,8 @@
 #include <algorithm>
 #include <sstream>
 
+#include <boost/tuple/tuple.hpp>
+
 #include <Standard_Failure.hxx>
 #include <Standard_ErrorHandler.hxx>
 
 using namespace std;
 using namespace SMESH::Controls;
 
-typedef map<const SMDS_MeshElement*, list<const SMDS_MeshNode*> >    TElemOfNodeListMap;
-typedef map<const SMDS_MeshElement*, list<const SMDS_MeshElement*> > TElemOfElemListMap;
-
 typedef SMDS_SetIterator< SMDS_pElement, TIDSortedElemSet::const_iterator> TSetIterator;
 
 //=======================================================================
@@ -120,6 +119,19 @@ SMESH_MeshEditor::SMESH_MeshEditor( SMESH_Mesh* theMesh )
 {
 }
 
+//================================================================================
+/*!
+ * \brief Clears myLastCreatedNodes and myLastCreatedElems
+ */
+//================================================================================
+
+void SMESH_MeshEditor::CrearLastCreated()
+{
+  myLastCreatedNodes.Clear();
+  myLastCreatedElems.Clear();
+}
+
+
 //=======================================================================
 /*!
  * \brief Add element
@@ -389,6 +401,44 @@ int SMESH_MeshEditor::Remove (const list< int >& theIDs,
   return removed;
 }
 
+//================================================================================
+/*!
+ * \brief Create 0D elements on all nodes of the given object except those
+ *        nodes on which a 0D element already exists.
+ *  \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
+ */
+//================================================================================
+
+void SMESH_MeshEditor::Create0DElementsOnAllNodes( const TIDSortedElemSet& elements,
+                                                   TIDSortedElemSet&       all0DElems )
+{
+  typedef SMDS_SetIterator<const SMDS_MeshElement*, TIDSortedElemSet::const_iterator> TSetIterator;
+  SMDS_ElemIteratorPtr elemIt;
+  if ( elements.empty() )
+    elemIt = GetMeshDS()->elementsIterator( SMDSAbs_Node );
+  else
+    elemIt = SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() ));
+
+  while ( elemIt->more() )
+  {
+    const SMDS_MeshElement* e = elemIt->next();
+    SMDS_ElemIteratorPtr nodeIt = e->nodesIterator();
+    while ( nodeIt->more() )
+    {
+      const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
+      SMDS_ElemIteratorPtr it0D = n->GetInverseElementIterator( SMDSAbs_0DElement );
+      if ( it0D->more() )
+        all0DElems.insert( it0D->next() );
+      else {
+        myLastCreatedElems.Append( GetMeshDS()->Add0DElement( n ));
+        all0DElems.insert( myLastCreatedElems.Last() );
+      }
+    }
+  }
+}
+
 //=======================================================================
 //function : FindShape
 //purpose  : Return an index of the shape theElem is on
@@ -1065,7 +1115,7 @@ 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 sould be orientated according to
+ * \param theFace - one of \a theFaces that sould be oriented according to
  *        \a theDirection and whose orientation defines orientation of other faces
  * \return number of reoriented faces.
  */
@@ -1093,7 +1143,7 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
 
   // Orient other faces
 
-  set< const SMDS_MeshElement* > startFaces;
+  set< const SMDS_MeshElement* > startFaces, visitedFaces;
   TIDSortedElemSet avoidSet;
   set< SMESH_TLink > checkedLinks;
   pair< set< SMESH_TLink >::iterator, bool > linkIt_isNew;
@@ -1102,16 +1152,26 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
     theFaces.erase( theFace );
   startFaces.insert( theFace );
 
+  int nodeInd1, nodeInd2;
+  const SMDS_MeshElement*           otherFace;
+  vector< const SMDS_MeshElement* > facesNearLink;
+  vector< std::pair< int, int > >   nodeIndsOfFace;
+
   set< const SMDS_MeshElement* >::iterator startFace = startFaces.begin();
-  while ( startFace != startFaces.end() )
+  while ( !startFaces.empty() )
   {
+    startFace = startFaces.begin();
     theFace = *startFace;
-    const int nbNodes = theFace->NbCornerNodes();
+    startFaces.erase( startFace );
+    if ( !visitedFaces.insert( theFace ).second )
+      continue;
 
     avoidSet.clear();
     avoidSet.insert(theFace);
 
     NLink link( theFace->GetNode( 0 ), 0 );
+
+    const int nbNodes = theFace->NbCornerNodes();
     for ( int i = 0; i < nbNodes; ++i ) // loop on links of theFace
     {
       link.second = theFace->GetNode(( i+1 ) % nbNodes );
@@ -1120,33 +1180,61 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
       {
         // link has already been checked and won't be encountered more
         // if the group (theFaces) is manifold
-        checkedLinks.erase( linkIt_isNew.first );
+        //checkedLinks.erase( linkIt_isNew.first );
       }
       else
       {
-        int nodeInd1, nodeInd2;
-        const SMDS_MeshElement* otherFace = FindFaceInSet( link.first, link.second,
-                                                           theFaces, avoidSet,
-                                                           & nodeInd1, & nodeInd2);
+        facesNearLink.clear();
+        nodeIndsOfFace.clear();
+        while (( otherFace = FindFaceInSet( link.first, link.second,
+                                            theFaces, avoidSet, &nodeInd1, &nodeInd2 )))
+          if ( otherFace != theFace)
+          {
+            facesNearLink.push_back( otherFace );
+            nodeIndsOfFace.push_back( make_pair( nodeInd1, nodeInd2 ));
+            avoidSet.insert( otherFace );
+          }
+        if ( facesNearLink.size() > 1 )
+        {
+          // NON-MANIFOLD mesh shell !
+          // select a face most co-directed with theFace,
+          // other faces won't be visited this time
+          gp_XYZ NF, NOF;
+          SMESH_Algo::FaceNormal( theFace, NF, /*normalized=*/false );
+          double proj, maxProj = -1;
+          for ( size_t i = 0; i < facesNearLink.size(); ++i ) {
+            SMESH_Algo::FaceNormal( facesNearLink[i], NOF, /*normalized=*/false );
+            if (( proj = Abs( NF * NOF )) > maxProj ) {
+              maxProj = proj;
+              otherFace = facesNearLink[i];
+              nodeInd1  = nodeIndsOfFace[i].first;
+              nodeInd2  = nodeIndsOfFace[i].second;
+            }
+          }
+          // 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] );
+        }
+        else if ( facesNearLink.size() == 1 )
+        {
+          otherFace = facesNearLink[0];
+          nodeInd1  = nodeIndsOfFace.back().first;
+          nodeInd2  = nodeIndsOfFace.back().second;
+        }
         if ( otherFace && otherFace != theFace)
         {
-          // link must be reversed in otherFace if orientation ot otherFace
+          // link must be reverse in otherFace if orientation ot otherFace
           // is same as that of theFace
           if ( abs(nodeInd2-nodeInd1) == 1 ? nodeInd2 > nodeInd1 : nodeInd1 > nodeInd2 )
           {
-            // cout << "Reorient " << otherFace->GetID() << " near theFace=" <<theFace->GetID()
-            //      << " \tlink( " << link.first->GetID() << " " << link.second->GetID() << endl;
             nbReori += Reorient( otherFace );
           }
           startFaces.insert( otherFace );
-          if ( theFaces.size() > 1 ) // leave 1 face to prevent finding not selected faces
-            theFaces.erase( otherFace );
         }
       }
-      std::swap( link.first, link.second );
+      std::swap( link.first, link.second ); // reverse the link
     }
-    startFaces.erase( startFace );
-    startFace = startFaces.begin();
   }
   return nbReori;
 }
@@ -1216,7 +1304,8 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet &                   theElems,
     if( !elem->IsQuadratic() ) {
 
       // split liner quadrangle
-
+      // for MaxElementLength2D functor we return minimum diagonal for splitting,
+      // because aBadRate1=2*len(diagonal 1-3); aBadRate2=2*len(diagonal 2-4)
       if ( aBadRate1 <= aBadRate2 ) {
         // tr1 + tr2 is better
         newElem1 = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
@@ -1350,7 +1439,8 @@ int SMESH_MeshEditor::BestSplit (const SMDS_MeshElement*              theQuad,
     SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
     SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
     aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
-
+    // for MaxElementLength2D functor we return minimum diagonal for splitting,
+    // because aBadRate1=2*len(diagonal 1-3); aBadRate2=2*len(diagonal 2-4)
     if (aBadRate1 <= aBadRate2) // tr1 + tr2 is better
       return 1; // diagonal 1-3
 
@@ -3062,7 +3152,7 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
     SMDS_FaceIteratorPtr fIt = aMesh->facesIterator();
     while ( fIt->more() ) {
       const SMDS_MeshElement* face = fIt->next();
-      theElems.insert( face );
+      theElems.insert( theElems.end(), face );
     }
   }
   // get all face ids theElems are on
@@ -4078,7 +4168,7 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
                                vecNewNodes[ 1 ]->second.back())) {
           myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
                                                    vecNewNodes[ 1 ]->second.back()));
-          srcElements.Append( myLastCreatedElems.Last() );
+          srcElements.Append( elem );
         }
       }
       else {
@@ -4088,7 +4178,7 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
           myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
                                                    vecNewNodes[ 1 ]->second.back(),
                                                    vecNewNodes[ 2 ]->second.back()));
-          srcElements.Append( myLastCreatedElems.Last() );
+          srcElements.Append( elem );
         }
       }
     }
@@ -4110,19 +4200,20 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
         int iNext = ( iNode + 1 == nbNodes ) ? 0 : iNode + 1;
         const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first;
         const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first;
-        // check if a link is free
+        // check if a link n1-n2 is free
         if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet )) {
           hasFreeLinks = true;
-          // make an edge and a ceiling for a new edge
-          if ( !aMesh->FindEdge( n1, n2 )) {
-            myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // free link edge
+          // make a new edge and a ceiling for a new edge
+          const SMDS_MeshElement* edge;
+          if ( ! ( edge = aMesh->FindEdge( n1, n2 ))) {
+            myLastCreatedElems.Append( edge = aMesh->AddEdge( n1, n2 )); // free link edge
             srcElements.Append( myLastCreatedElems.Last() );
           }
           n1 = vecNewNodes[ iNode ]->second.back();
           n2 = vecNewNodes[ iNext ]->second.back();
           if ( !aMesh->FindEdge( n1, n2 )) {
-            myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // ceiling edge
-            srcElements.Append( myLastCreatedElems.Last() );
+            myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // new edge ceiling
+            srcElements.Append( edge );
           }
         }
       }
@@ -4144,14 +4235,14 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
           // find medium node
           if ( !aMesh->FindEdge( n1, n2, n3 )) {
             myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // free link edge
-            srcElements.Append( myLastCreatedElems.Last() );
+            srcElements.Append( elem );
           }
           n1 = vecNewNodes[ iNode ]->second.back();
           n2 = vecNewNodes[ iNext ]->second.back();
           n3 = vecNewNodes[ iNode+nbn ]->second.back();
           if ( !aMesh->FindEdge( n1, n2, n3 )) {
             myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // ceiling edge
-            srcElements.Append( myLastCreatedElems.Last() );
+            srcElements.Append( elem );
           }
         }
       }
@@ -4407,7 +4498,7 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
       }
 
       while ( srcElements.Length() < myLastCreatedElems.Length() )
-        srcElements.Append( myLastCreatedElems.Last() );
+        srcElements.Append( elem );
     }
   } // loop on swept elements
 }
@@ -5082,10 +5173,16 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
     if ( BRep_Tool::Degenerated( aTrackEdge ) )
       return EXTR_BAD_PATH_SHAPE;
     TopExp::Vertices( aTrackEdge, aV1, aV2 );
-    aItN = theTrack->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
-    const SMDS_MeshNode* aN1 = aItN->next();
-    aItN = theTrack->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
-    const SMDS_MeshNode* aN2 = aItN->next();
+    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();
+    }
     // starting node must be aN1 or aN2
     if ( !( aN1 == theN1 || aN2 == theN1 ) )
       return EXTR_BAD_STARTING_NODE;
@@ -5115,7 +5212,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
       }
     }
     list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
-    int startNid = theN1->GetID();
+    TopoDS_Vertex aVprev;
     TColStd_MapOfInteger UsedNums;
     int NbEdges = Edges.Length();
     int i = 1;
@@ -5129,17 +5226,37 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
         SMESH_subMesh* locTrack = *itLSM;
         SMESHDS_SubMesh* locMeshDS = locTrack->GetSubMeshDS();
         TopExp::Vertices( aTrackEdge, aV1, aV2 );
-        aItN = locTrack->GetFather()->GetSubMesh(aV1)->GetSubMeshDS()->GetNodes();
-        const SMDS_MeshNode* aN1 = aItN->next();
-        aItN = locTrack->GetFather()->GetSubMesh(aV2)->GetSubMeshDS()->GetNodes();
-        const SMDS_MeshNode* aN2 = aItN->next();
-        // starting node must be aN1 or aN2
-        if ( !( aN1->GetID() == startNid || aN2->GetID() == startNid ) ) continue;
+        bool aN1isOK = false, aN2isOK = false;
+        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();
+          }
+          // starting node must be aN1 or aN2
+          aN1isOK = ( aN1 && aN1 == theN1 );
+          aN2isOK = ( aN2 && aN2 == theN1 );
+        }
+        else {
+          // we have specified ending vertex of the previous edge on the previous iteration
+          // and we have just to check that it corresponds to any vertex in current segment
+          aN1isOK = aVprev.IsSame( aV1 );
+          aN2isOK = aVprev.IsSame( aV2 );
+        }
+        if ( !aN1isOK && !aN2isOK ) continue;
         // 2. Collect parameters on the track edge
         aPrms.clear();
         aItN = locMeshDS->GetNodes();
         while ( aItN->more() ) {
-          const SMDS_MeshNode* pNode = aItN->next();
+          const SMDS_MeshNode*     pNode = aItN->next();
           const SMDS_EdgePosition* pEPos =
             static_cast<const SMDS_EdgePosition*>( pNode->GetPosition() );
           double aT = pEPos->GetUParameter();
@@ -5147,12 +5264,12 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
         }
         list<SMESH_MeshEditor_PathPoint> LPP;
         //Extrusion_Error err =
-        MakeEdgePathPoints(aPrms, aTrackEdge,(aN1->GetID()==startNid), LPP);
+        MakeEdgePathPoints(aPrms, aTrackEdge, aN1isOK, LPP);
         LLPPs.push_back(LPP);
         UsedNums.Add(k);
         // update startN for search following egde
-        if( aN1->GetID() == startNid ) startNid = aN2->GetID();
-        else startNid = aN1->GetID();
+        if ( aN1isOK ) aVprev = aV2;
+        else           aVprev = aV1;
         break;
       }
     }
@@ -5171,8 +5288,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet &   theElements,
       SMESH_MeshEditor_PathPoint PP2 = currList.front();
       gp_Dir D1 = PP1.Tangent();
       gp_Dir D2 = PP2.Tangent();
-      gp_Dir Dnew( gp_Vec( (D1.X()+D2.X())/2, (D1.Y()+D2.Y())/2,
-                           (D1.Z()+D2.Z())/2 ) );
+      gp_Dir Dnew( ( D1.XYZ() + D2.XYZ() ) / 2 );
       PP1.SetTangent(Dnew);
       fullList.push_back(PP1);
       itPP++;
@@ -5885,13 +6001,15 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
 
   // Sort existing groups by types and collect their names
 
-  // to store an old group and a generated new one
-  typedef pair< SMESHDS_GroupBase*, SMESHDS_Group* > TOldNewGroup;
+  // to store an old group and a generated new ones
+  using boost::tuple;
+  using boost::make_tuple;
+  typedef tuple< SMESHDS_GroupBase*, SMESHDS_Group*, SMESHDS_Group* > TOldNewGroup;
   vector< list< TOldNewGroup > > groupsByType( SMDSAbs_NbElementTypes );
   vector< TOldNewGroup* > orderedOldNewGroups; // in order of old groups
   // group names
   set< string > groupNames;
-  
+
   SMESH_Mesh::GroupIteratorPtr groupIt = GetMesh()->GetGroups();
   if ( !groupIt->more() ) return newGroupIDs;
 
@@ -5902,12 +6020,13 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
     if ( !group ) continue;
     SMESHDS_GroupBase* groupDS = group->GetGroupDS();
     if ( !groupDS || groupDS->IsEmpty() ) continue;
-    groupNames.insert( group->GetName() );
+    groupNames.insert    ( group->GetName() );
     groupDS->SetStoreName( group->GetName() );
-    SMESHDS_Group* newGroup = new SMESHDS_Group( newGroupID++, mesh->GetMeshDS(),
-                                                 groupDS->GetType() );
-    groupsByType[ groupDS->GetType() ].push_back( make_pair( groupDS, newGroup ));
-    orderedOldNewGroups.push_back( & groupsByType[ groupDS->GetType() ].back() );
+    const SMDSAbs_ElementType type = groupDS->GetType();
+    SMESHDS_Group* newGroup    = new SMESHDS_Group( newGroupID++, mesh->GetMeshDS(), type );
+    SMESHDS_Group* newTopGroup = new SMESHDS_Group( newGroupID++, mesh->GetMeshDS(), type );
+    groupsByType[ type ].push_back( make_tuple( groupDS, newGroup, newTopGroup ));
+    orderedOldNewGroups.push_back( & groupsByType[ type ].back() );
   }
 
   // Loop on nodes and elements to add them in new groups
@@ -5917,7 +6036,7 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
     const SMESH_SequenceOfElemPtr& gens  = isNodes ? nodeGens : elemGens;
     const SMESH_SequenceOfElemPtr& elems = isNodes ? myLastCreatedNodes : myLastCreatedElems;
     if ( gens.Length() != elems.Length() )
-      throw SALOME_Exception(LOCALIZED("invalid args"));
+      throw SALOME_Exception("SMESH_MeshEditor::generateGroups(): invalid args");
 
     // loop on created elements
     for (int iElem = 1; iElem <= elems.Length(); ++iElem )
@@ -5933,7 +6052,7 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
           ++iElem; // skip all elements made by sourceElem
         continue;
       }
-      // collect all elements made by sourceElem
+      // collect all elements made by the iElem-th sourceElem
       list< const SMDS_MeshElement* > resultElems;
       if ( const SMDS_MeshElement* resElem = elems( iElem ))
         if ( resElem != sourceElem )
@@ -5943,18 +6062,44 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
           if ( resElem != sourceElem )
             resultElems.push_back( resElem );
 
-      // add resultElems to groups made by ones the sourceElem belongs to
+      // there must be a top element
+      const SMDS_MeshElement* topElem = 0;
+      if ( isNodes )
+      {
+        topElem = resultElems.back();
+        resultElems.pop_back();
+      }
+      else
+      {
+        list< const SMDS_MeshElement* >::reverse_iterator resElemIt = resultElems.rbegin();
+        for ( ; resElemIt != resultElems.rend() ; ++resElemIt )
+          if ( (*resElemIt)->GetType() == sourceElem->GetType() )
+          {
+            topElem = *resElemIt;
+            resultElems.erase( --(resElemIt.base()) ); // erase *resElemIt
+            break;
+          }
+      }
+
+      // add resultElems to groups originted from ones the sourceElem belongs to
       list< TOldNewGroup >::iterator gOldNew, gLast = groupsOldNew.end();
       for ( gOldNew = groupsOldNew.begin(); gOldNew != gLast; ++gOldNew )
       {
-        SMESHDS_GroupBase* oldGroup = gOldNew->first;
+        SMESHDS_GroupBase* oldGroup = gOldNew->get<0>();
         if ( oldGroup->Contains( sourceElem )) // sourceElem is in oldGroup
         {
           // fill in a new group
-          SMDS_MeshGroup & newGroup = gOldNew->second->SMDSGroup();
+          SMDS_MeshGroup & newGroup = gOldNew->get<1>()->SMDSGroup();
           list< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt;
           for ( resElemIt = resultElems.begin(); resElemIt != resLast; ++resElemIt )
             newGroup.Add( *resElemIt );
+
+          // fill a "top" group
+          if ( topElem )
+          {
+            SMDS_MeshGroup & newTopGroup = gOldNew->get<2>()->SMDSGroup();
+            newTopGroup.Add( topElem );
+          }
         }
       }
     } // loop on created elements
@@ -5962,35 +6107,54 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
 
   // Create new SMESH_Groups from SMESHDS_Groups and remove empty SMESHDS_Groups
 
+  list<int> topGrouIds;
   for ( size_t i = 0; i < orderedOldNewGroups.size(); ++i )
   {
-    SMESHDS_GroupBase* oldGroupDS = orderedOldNewGroups[i]->first;
-    SMESHDS_Group*     newGroupDS = orderedOldNewGroups[i]->second;
-    if ( newGroupDS->IsEmpty() )
-    {
-      mesh->GetMeshDS()->RemoveGroup( newGroupDS );
-    }
-    else
-    {
-      // make a name
-      string name = oldGroupDS->GetStoreName();
-      if ( !targetMesh ) {
-        name += "_";
-        name += postfix;
-        int nb = 1;
-        while ( !groupNames.insert( name ).second ) // name exists
-          name = SMESH_Comment( oldGroupDS->GetStoreName() ) << "_" << postfix << "_" << nb++;
+    SMESHDS_GroupBase* oldGroupDS =   orderedOldNewGroups[i]->get<0>();
+    SMESHDS_Group*   newGroups[2] = { orderedOldNewGroups[i]->get<1>(),
+                                      orderedOldNewGroups[i]->get<2>() };
+    const int nbNewGroups = !newGroups[0]->IsEmpty() + !newGroups[1]->IsEmpty();
+    for ( int is2nd = 0; is2nd < 2; ++is2nd )
+    {
+      SMESHDS_Group* newGroupDS = newGroups[ is2nd ];
+      if ( newGroupDS->IsEmpty() )
+      {
+        mesh->GetMeshDS()->RemoveGroup( newGroupDS );
+      }
+      else
+      {
+        // set group type
+        newGroupDS->SetType( newGroupDS->GetElements()->next()->GetType() );
+
+        // make a name
+        const bool isTop = ( nbNewGroups == 2 &&
+                             newGroupDS->GetType() == oldGroupDS->GetType() &&
+                             is2nd );
+
+        string name = oldGroupDS->GetStoreName();
+        if ( !targetMesh ) {
+          string suffix = ( isTop ? "top": postfix.c_str() );
+          name += "_";
+          name += suffix;
+          int nb = 1;
+          while ( !groupNames.insert( name ).second ) // name exists
+            name = SMESH_Comment( oldGroupDS->GetStoreName() ) << "_" << suffix << "_" << nb++;
+        }
+        else if ( isTop ) {
+          name += "_top";
+        }
+        newGroupDS->SetStoreName( name.c_str() );
+
+        // make a SMESH_Groups
+        mesh->AddGroup( newGroupDS );
+        if ( isTop )
+          topGrouIds.push_back( newGroupDS->GetID() );
+        else
+          newGroupIDs->push_back( newGroupDS->GetID() );
       }
-      newGroupDS->SetStoreName( name.c_str() );
-
-      // make a SMESH_Groups
-      mesh->AddGroup( newGroupDS );
-      newGroupIDs->push_back( newGroupDS->GetID() );
-
-      // set group type
-      newGroupDS->SetType( newGroupDS->GetElements()->next()->GetType() );
     }
   }
+  newGroupIDs->splice( newGroupIDs->end(), topGrouIds );
 
   return newGroupIDs;
 }
@@ -6103,7 +6267,7 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
         }
         else if ( tree->NbNodes() ) // put a tree to the treeMap
         {
-          const Bnd_B3d& box = tree->getBox();
+          const Bnd_B3d& box = *tree->getBox();
           double sqDist = thePnt.SquareDistance( 0.5 * ( box.CornerMin() + box.CornerMax() ));
           pair<TDistTreeMap::iterator,bool> it_in = treeMap.insert( make_pair( sqDist, tree ));
           if ( !it_in.second ) // not unique distance to box center
@@ -6115,7 +6279,7 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
       TDistTreeMap::iterator sqDist_tree = treeMap.begin();
       if ( treeMap.size() > 5 ) {
         SMESH_OctreeNode* closestTree = sqDist_tree->second;
-        const Bnd_B3d& box = closestTree->getBox();
+        const Bnd_B3d& box = *closestTree->getBox();
         double limit = sqrt( sqDist_tree->first ) + sqrt ( box.SquareExtent() );
         sqLimit = limit * limit;
       }
@@ -6200,7 +6364,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
 
   protected:
     ElementBndBoxTree():_size(0) {}
-    SMESH_Octree* allocateOctreeChild() const { return new ElementBndBoxTree; }
+    SMESH_Octree* newChild() const { return new ElementBndBoxTree; }
     void          buildChildrenData();
     Bnd_B3d*      buildRootBox();
   private:
@@ -6222,7 +6386,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
   //================================================================================
 
   ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, SMDS_ElemIteratorPtr theElemIt, double tolerance)
-    :SMESH_Octree( new SMESH_Octree::Limit( MaxLevel, /*minSize=*/0. ))
+    :SMESH_Octree( new SMESH_TreeLimit( MaxLevel, /*minSize=*/0. ))
   {
     int nbElems = mesh.GetMeshInfo().NbElements( elemType );
     _elements.reserve( nbElems );
@@ -6273,7 +6437,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
     {
       for (int j = 0; j < 8; j++)
       {
-        if ( !_elements[i]->IsOut( myChildren[j]->getBox() ))
+        if ( !_elements[i]->IsOut( *myChildren[j]->getBox() ))
         {
           _elements[i]->_refCount++;
           ((ElementBndBoxTree*)myChildren[j])->_elements.push_back( _elements[i]);
@@ -6304,7 +6468,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
   void ElementBndBoxTree::getElementsNearPoint( const gp_Pnt&     point,
                                                 TIDSortedElemSet& foundElems)
   {
-    if ( getBox().IsOut( point.XYZ() ))
+    if ( getBox()->IsOut( point.XYZ() ))
       return;
 
     if ( isLeaf() )
@@ -6329,7 +6493,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
   void ElementBndBoxTree::getElementsNearLine( const gp_Ax1&     line,
                                                TIDSortedElemSet& foundElems)
   {
-    if ( getBox().IsOut( line ))
+    if ( getBox()->IsOut( line ))
       return;
 
     if ( isLeaf() )
@@ -6355,7 +6519,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
                                                 const double      radius,
                                                 TIDSortedElemSet& foundElems)
   {
-    if ( getBox().IsOut( center, radius ))
+    if ( getBox()->IsOut( center, radius ))
       return;
 
     if ( isLeaf() )
@@ -6657,7 +6821,7 @@ void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerF
  * \brief Find elements of given type where the given point is IN or ON.
  *        Returns nb of found elements and elements them-selves.
  *
- * 'ALL' type means elements of any type excluding nodes, balls and 0D elements 
+ * 'ALL' type means elements of any type excluding nodes, balls and 0D elements
  */
 //=======================================================================
 
@@ -6734,13 +6898,13 @@ SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt&       point,
     }
     TIDSortedElemSet suspectElems;
     _ebbTree->getElementsNearPoint( point, suspectElems );
-    
+
     if ( suspectElems.empty() && _ebbTree->maxSize() > 0 )
     {
-      gp_Pnt boxCenter = 0.5 * ( _ebbTree->getBox().CornerMin() +
-                                 _ebbTree->getBox().CornerMax() );
+      gp_Pnt boxCenter = 0.5 * ( _ebbTree->getBox()->CornerMin() +
+                                 _ebbTree->getBox()->CornerMax() );
       double radius;
-      if ( _ebbTree->getBox().IsOut( point.XYZ() ))
+      if ( _ebbTree->getBox()->IsOut( point.XYZ() ))
         radius = point.Distance( boxCenter ) - 0.5 * _ebbTree->maxSize();
       else
         radius = _ebbTree->maxSize() / pow( 2., _ebbTree->getHeight()) / 2;
@@ -7242,7 +7406,7 @@ namespace
                       POS_ALL = POS_LEFT | POS_RIGHT | POS_VERTEX };
   struct PointPos
   {
-    PositionName _name; 
+    PositionName _name;
     int          _index; // index of vertex or segment
 
     PointPos( PositionName n, int i=-1 ): _name(n), _index(i) {}
@@ -7506,6 +7670,12 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
         //MESSAGE("  node to remove " << nToRemove->GetID());
         rmNodeIds.push_back( nToRemove->GetID() );
         AddToSameGroups( nToKeep, nToRemove, aMesh );
+        // set _alwaysComputed to a sub-mesh of VERTEX to enable mesh computing
+        // after MergeNodes() w/o creating node in place of merged ones.
+        const SMDS_PositionPtr& pos = nToRemove->GetPosition();
+        if ( pos && pos->GetTypeOfPosition() == SMDS_TOP_VERTEX )
+          if ( SMESH_subMesh* sm = myMesh->GetSubMeshContaining( nToRemove->getshapeId() ))
+            sm->SetIsAlwaysComputed( true );
       }
 
       SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
@@ -8156,32 +8326,29 @@ private:
 //purpose  : Return list of group of elements built on the same nodes.
 //           Search among theElements or in the whole mesh if theElements is empty
 //=======================================================================
-void SMESH_MeshEditor::FindEqualElements(set<const SMDS_MeshElement*> & theElements,
-                                         TListOfListOfElementsID &      theGroupsOfElementsID)
+
+void SMESH_MeshEditor::FindEqualElements(TIDSortedElemSet &        theElements,
+                                         TListOfListOfElementsID & theGroupsOfElementsID)
 {
   myLastCreatedElems.Clear();
   myLastCreatedNodes.Clear();
 
-  typedef set<const SMDS_MeshElement*> TElemsSet;
   typedef map< SortableElement, int > TMapOfNodeSet;
   typedef list<int> TGroupOfElems;
 
-  TElemsSet elems;
   if ( theElements.empty() )
   { // get all elements in the mesh
     SMDS_ElemIteratorPtr eIt = GetMeshDS()->elementsIterator();
     while ( eIt->more() )
-      elems.insert( elems.end(), eIt->next());
+      theElements.insert( theElements.end(), eIt->next());
   }
-  else
-    elems = theElements;
 
   vector< TGroupOfElems > arrayOfGroups;
   TGroupOfElems groupOfElems;
   TMapOfNodeSet mapOfNodeSet;
 
-  TElemsSet::iterator elemIt = elems.begin();
-  for ( int i = 0, j=0; elemIt != elems.end(); ++elemIt, ++j ) {
+  TIDSortedElemSet::iterator elemIt = theElements.begin();
+  for ( int i = 0, j=0; elemIt != theElements.end(); ++elemIt, ++j ) {
     const SMDS_MeshElement* curElem = *elemIt;
     SortableElement SE(curElem);
     int ind = -1;
@@ -8254,8 +8421,8 @@ void SMESH_MeshEditor::MergeElements(TListOfListOfElementsID & theGroupsOfElemen
 
 void SMESH_MeshEditor::MergeEqualElements()
 {
-  set<const SMDS_MeshElement*> aMeshElements; /* empty input -
-                                                 to merge equal elements in the whole mesh */
+  TIDSortedElemSet aMeshElements; /* empty input ==
+                                     to merge equal elements in the whole mesh */
   TListOfListOfElementsID aGroupsOfElementsID;
   FindEqualElements(aMeshElements, aGroupsOfElementsID);
   MergeElements(aGroupsOfElementsID);
@@ -9440,14 +9607,26 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
   {
     nbElem++;
     const SMDS_MeshElement* elem = ElemItr->next();
-    if( !elem || elem->IsQuadratic() ) continue;
+    if( !elem ) continue;
 
+    const SMDSAbs_EntityType aGeomType = elem->GetEntityType();
+    if ( elem->IsQuadratic() )
+    {
+      bool alreadyOK;
+      switch ( aGeomType ) {
+      case SMDSEntity_Quad_Quadrangle:
+      case SMDSEntity_Quad_Hexa:         alreadyOK = !theHelper.GetIsBiQuadratic(); break;
+      case SMDSEntity_BiQuad_Quadrangle:
+      case SMDSEntity_TriQuad_Hexa:      alreadyOK = theHelper.GetIsBiQuadratic(); break;
+      default:                           alreadyOK = true;
+      }
+      if ( alreadyOK ) continue;
+    }
     // get elem data needed to re-create it
     //
-    const int id                        = elem->GetID();
-    const int nbNodes                   = elem->NbNodes();
-    const SMDSAbs_ElementType aType     = elem->GetType();
-    const SMDSAbs_EntityType  aGeomType = elem->GetEntityType();
+    const int id                     = elem->GetID();
+    const int nbNodes                = elem->NbCornerNodes();
+    const SMDSAbs_ElementType aType  = elem->GetType();
     nodes.assign(elem->begin_nodes(), elem->end_nodes());
     if ( aGeomType == SMDSEntity_Polyhedra )
       nbNodeInFaces = static_cast<const SMDS_VtkVolume* >( elem )->GetQuantities();
@@ -9496,6 +9675,8 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
           NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], nodes[5], id, theForce3d);
           break;
         case SMDSEntity_Hexa:
+        case SMDSEntity_Quad_Hexa:
+        case SMDSEntity_TriQuad_Hexa:
           NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
                                         nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
           break;
@@ -9514,18 +9695,20 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
   }
   return nbElem;
 }
-
 //=======================================================================
 //function : ConvertToQuadratic
 //purpose  :
 //=======================================================================
 
-void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
+void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theToBiQuad)
 {
   SMESHDS_Mesh* meshDS = GetMeshDS();
 
   SMESH_MesherHelper aHelper(*myMesh);
+
   aHelper.SetIsQuadratic( true );
+  aHelper.SetIsBiQuadratic( theToBiQuad );
+  aHelper.SetElementsOnShape(true);
 
   int nbCheckedElems = 0;
   if ( myMesh->HasShapeToMesh() )
@@ -9545,6 +9728,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
   int totalNbElems = meshDS->NbEdges() + meshDS->NbFaces() + meshDS->NbVolumes();
   if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes
   {
+    aHelper.SetElementsOnShape(false);
     SMESHDS_SubMesh *smDS = 0;
     SMDS_EdgeIteratorPtr aEdgeItr = meshDS->edgesIterator();
     while(aEdgeItr->more())
@@ -9567,10 +9751,14 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
     while(aFaceItr->more())
     {
       const SMDS_MeshFace* face = aFaceItr->next();
-      if(!face || face->IsQuadratic() ) continue;
+      if ( !face ) continue;
+      
+      const SMDSAbs_EntityType type = face->GetEntityType();
+      if (( theToBiQuad  && type == SMDSEntity_BiQuad_Quadrangle ) ||
+          ( !theToBiQuad && type == SMDSEntity_Quad_Quadrangle ))
+        continue;
 
       const int id = face->GetID();
-      const SMDSAbs_EntityType type = face->GetEntityType();
       vector<const SMDS_MeshNode *> nodes ( face->begin_nodes(), face->end_nodes());
 
       meshDS->RemoveFreeElement(face, smDS, /*fromGroups=*/false);
@@ -9596,8 +9784,12 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
       const SMDS_MeshVolume* volume = aVolumeItr->next();
       if(!volume || volume->IsQuadratic() ) continue;
 
-      const int id = volume->GetID();
       const SMDSAbs_EntityType type = volume->GetEntityType();
+      if (( theToBiQuad  && type == SMDSEntity_TriQuad_Hexa ) ||
+          ( !theToBiQuad && type == SMDSEntity_Quad_Hexa ))
+        continue;
+
+      const int id = volume->GetID();
       vector<const SMDS_MeshNode *> nodes (volume->begin_nodes(), volume->end_nodes());
       if ( type == SMDSEntity_Polyhedra )
         nbNodeInFaces = static_cast<const SMDS_VtkVolume* >(volume)->GetQuantities();
@@ -9613,6 +9805,8 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
         NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d );
         break;
       case SMDSEntity_Hexa:
+      case SMDSEntity_Quad_Hexa:
+      case SMDSEntity_TriQuad_Hexa:
         NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
                                       nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
         break;
@@ -9648,7 +9842,8 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
 //================================================================================
 
 void SMESH_MeshEditor::ConvertToQuadratic(const bool        theForce3d,
-                                          TIDSortedElemSet& theElements)
+                                          TIDSortedElemSet& theElements,
+                                          const bool        theToBiQuad)
 {
   if ( theElements.empty() ) return;
 
@@ -9675,8 +9870,19 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool        theForce3d,
       const SMDS_MeshElement* e = invIt->next();
       if ( e->IsQuadratic() )
       {
-        quadAdjacentElems[ e->GetType() ].insert( e );
-        continue;
+        bool alreadyOK;
+        switch ( e->GetEntityType() ) {
+        case SMDSEntity_Quad_Quadrangle:
+        case SMDSEntity_Quad_Hexa:         alreadyOK = !theToBiQuad; break;
+        case SMDSEntity_BiQuad_Quadrangle:
+        case SMDSEntity_TriQuad_Hexa:      alreadyOK = theToBiQuad; break;
+        default:                           alreadyOK = true;
+        }
+        if ( alreadyOK )
+        {
+          quadAdjacentElems[ e->GetType() ].insert( e );
+          continue;
+        }
       }
       if ( e->GetType() >= elemType )
       {
@@ -9698,6 +9904,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool        theForce3d,
 
   SMESH_MesherHelper helper(*myMesh);
   helper.SetIsQuadratic( true );
+  helper.SetIsBiQuadratic( theToBiQuad );
 
   // add links of quadratic adjacent elements to the helper
 
@@ -9720,18 +9927,32 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool        theForce3d,
       helper.AddTLinks( static_cast< const SMDS_MeshVolume*> (*eIt) );
     }
 
-  // make quadratic elements instead of linear ones
+  // make quadratic (or bi-tri-quadratic) elements instead of linear ones
 
-  SMESHDS_Mesh* meshDS = GetMeshDS();
+  SMESHDS_Mesh*  meshDS = GetMeshDS();
   SMESHDS_SubMesh* smDS = 0;
   for ( eIt = theElements.begin(); eIt != theElements.end(); ++eIt )
   {
     const SMDS_MeshElement* elem = *eIt;
-    if( elem->IsQuadratic() || elem->NbNodes() < 2 || elem->IsPoly() )
+    if( elem->NbNodes() < 2 || elem->IsPoly() )
       continue;
 
-    const int id                   = elem->GetID();
+    if ( elem->IsQuadratic() )
+    {
+      bool alreadyOK;
+      switch ( elem->GetEntityType() ) {
+      case SMDSEntity_Quad_Quadrangle:
+      case SMDSEntity_Quad_Hexa:         alreadyOK = !theToBiQuad; break;
+      case SMDSEntity_BiQuad_Quadrangle:
+      case SMDSEntity_TriQuad_Hexa:      alreadyOK = theToBiQuad; break;
+      default:                           alreadyOK = true;
+      }
+      if ( alreadyOK ) continue;
+    }
+
     const SMDSAbs_ElementType type = elem->GetType();
+    const int                   id = elem->GetID();
+    const int              nbNodes = elem->NbCornerNodes();
     vector<const SMDS_MeshNode *> nodes ( elem->begin_nodes(), elem->end_nodes());
 
     if ( !smDS || !smDS->Contains( elem ))
@@ -9739,7 +9960,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool        theForce3d,
     meshDS->RemoveFreeElement(elem, smDS, /*fromGroups=*/false);
 
     SMDS_MeshElement * newElem = 0;
-    switch( nodes.size() )
+    switch( nbNodes )
     {
     case 4: // cases for most frequently used element types go first (for optimization)
       if ( type == SMDSAbs_Volume )
@@ -11044,11 +11265,14 @@ double SMESH_MeshEditor::OrientedAngle(const gp_Pnt& p0, const gp_Pnt& p1, const
 
 /*!
  * \brief Double nodes on shared faces between groups of volumes and create flat elements on demand.
- * The list of groups must describe a partition of the mesh volumes.
- * The nodes of the internal faces at the boundaries of the groups are doubled.
- * In option, the internal faces are replaced by flat elements.
- * Triangles are transformed in prisms, and quadrangles in hexahedrons.
- * The flat elements are stored in groups of volumes.
+ *  The list of groups must contain at least two groups. The groups have to be disjoint: no common element into two different groups.
+ * The nodes of the internal faces at the boundaries of the groups are doubled. Optionally, the internal faces are replaced by flat elements.
+ * Triangles are transformed into prisms, and quadrangles into hexahedrons.
+ * The flat elements are stored in groups of volumes. These groups are named according to the position of the group in the list:
+ * the group j_n_p is the group of the flat elements that are built between the group #n and the group #p in the list.
+ * If there is no shared faces between the group #n and the group #p in the list, the group j_n_p is not created.
+ * All the flat elements are gathered into the group named "joints3D" (or "joints2D" in 2D situation).
+ * The flat element of the multiple junctions between the simple junction are stored in a group named "jointsMultiples".
  * @param theElems - list of groups of volumes, where a group of volume is a set of
  * SMDS_MeshElements sorted by Id.
  * @param createJointElems - if TRUE, create the elements
@@ -11082,6 +11306,32 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   std::set<int> emptySet;
   emptyMap.clear();
 
+  MESSAGE(".. Number of domains :"<<theElems.size());
+
+  // Check if the domains do not share an element
+  for (int idom = 0; idom < theElems.size()-1; idom++)
+    {
+//       MESSAGE("... Check of domain #" << idom);
+      const TIDSortedElemSet& domain = theElems[idom];
+      TIDSortedElemSet::const_iterator elemItr = domain.begin();
+      for (; elemItr != domain.end(); ++elemItr)
+        {
+          SMDS_MeshElement* anElem = (SMDS_MeshElement*) *elemItr;
+          int idombisdeb = idom + 1 ;
+          for (int idombis = idombisdeb; idombis < theElems.size(); idombis++) // check if the element belongs to a domain further in the list
+          {
+            const TIDSortedElemSet& domainbis = theElems[idombis];
+            if ( domainbis.count(anElem) )
+            {
+              MESSAGE(".... Domain #" << idom);
+              MESSAGE(".... Domain #" << idombis);
+              throw SALOME_Exception("The domains are not disjoint.");
+              return false ;
+            }
+          }
+        }
+    }
+
   for (int idom = 0; idom < theElems.size(); idom++)
     {
 
@@ -11090,7 +11340,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
       //     and corresponding volume of this domain, for each shared face.
       //     a volume has a face shared by 2 domains if it has a neighbor which is not in his domain.
 
-      //MESSAGE("Domain " << idom);
+      MESSAGE("... Neighbors of domain #" << idom);
       const TIDSortedElemSet& domain = theElems[idom];
       TIDSortedElemSet::const_iterator elemItr = domain.begin();
       for (; elemItr != domain.end(); ++elemItr)
@@ -11110,15 +11360,25 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
               const SMDS_MeshElement* elem = meshDS->FindElement(smdsId);
               if (! domain.count(elem)) // neighbor is in another domain : face is shared
                 {
-                  DownIdType face(downIds[n], downTypes[n]);
-                  if (!faceDomains.count(face))
-                    faceDomains[face] = emptyMap; // create an empty entry for face
-                  if (!faceDomains[face].count(idom))
-                    {
-                      faceDomains[face][idom] = vtkId; // volume associated to face in this domain
-                      celldom[vtkId] = idom;
-                      //MESSAGE("       cell with a border " << vtkId << " domain " << idom);
-                    }
+                  bool ok = false ;
+                  for (int idombis = 0; idombis < theElems.size(); idombis++) // check if the neighbor belongs to another domain of the list
+                  {
+                    // MESSAGE("Domain " << idombis);
+                    const TIDSortedElemSet& domainbis = theElems[idombis];
+                    if ( domainbis.count(elem)) ok = true ; // neighbor is in a correct domain : face is kept
+                  }
+                  if ( ok ) // the characteristics of the face is stored
+                  {
+                    DownIdType face(downIds[n], downTypes[n]);
+                    if (!faceDomains.count(face))
+                      faceDomains[face] = emptyMap; // create an empty entry for face
+                    if (!faceDomains[face].count(idom))
+                      {
+                        faceDomains[face][idom] = vtkId; // volume associated to face in this domain
+                        celldom[vtkId] = idom;
+                        //MESSAGE("       cell with a border " << vtkId << " domain " << idom);
+                      }
+                  }
                 }
             }
         }
@@ -11189,6 +11449,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   std::map<int, std::vector<int> > mutipleNodes; // nodes multi domains with domain order
   std::map<int, std::vector<int> > mutipleNodesToFace; // nodes multi domains with domain order to transform in Face (junction between 3 or more 2D domains)
 
+  MESSAGE(".. Duplication of the nodes");
   for (int idomain = 0; idomain < theElems.size(); idomain++)
     {
       itface = faceDomains.begin();
@@ -11251,6 +11512,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
         }
     }
 
+  MESSAGE(".. Creation of elements");
   for (int idomain = 0; idomain < theElems.size(); idomain++)
     {
       itface = faceDomains.begin();
@@ -11380,6 +11642,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   std::map<int, std::map<long,int> > nodeQuadDomains;
   std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
 
+  MESSAGE(".. Creation of elements: simple junction");
   if (createJointElems)
     {
       int idg;
@@ -11430,6 +11693,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   //     iterate on mutipleNodesToFace
   //     iterate on edgesMultiDomains
 
+  MESSAGE(".. Creation of elements: multiple junction");
   if (createJointElems)
     {
       // --- iterate on mutipleNodesToFace
@@ -11502,6 +11766,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   faceOrEdgeDom.clear();
   feDom.clear();
 
+  MESSAGE(".. Modification of elements");
   for (int idomain = 0; idomain < theElems.size(); idomain++)
     {
       std::map<int, std::map<int, int> >::const_iterator itnod = nodeDomains.begin();
@@ -11779,8 +12044,11 @@ void SMESH_MeshEditor::CreateHoleSkin(double radius,
       groupDS = group->GetGroupDS();
       if ( !groupDS || groupDS->IsEmpty() ) continue;
       std::string grpName = group->GetName();
+      //MESSAGE("grpName=" << grpName);
       if (grpName == groupName)
         break;
+      else
+        groupDS = 0;
     }
 
   bool isNodeGroup = false;
@@ -11794,6 +12062,7 @@ void SMESH_MeshEditor::CreateHoleSkin(double radius,
 
   if (nodesCoords.size() > 0)
     isNodeCoords = true; // a list o nodes given by their coordinates
+  //MESSAGE("---" << isNodeGroup << " " << isNodeCoords);
 
   // --- define groups to build