Salome HOME
PR: debug generate flat elements
[modules/smesh.git] / src / SMESH / SMESH_MeshEditor.cxx
index 5579ca090397014aa5c8ebe76471ea91bf1fa51e..49f9c1b3a8776abf98730188553103d7e1ca84fb 100644 (file)
@@ -94,6 +94,8 @@
 #include <set>
 #include <numeric>
 #include <limits>
+#include <algorithm>
+#include <sstream>
 
 #define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
 
@@ -132,114 +134,125 @@ SMESH_MeshEditor::AddElement(const vector<const SMDS_MeshNode*> & node,
   int nbnode = node.size();
   SMESHDS_Mesh* mesh = GetMeshDS();
   switch ( type ) {
-  case SMDSAbs_0DElement:
-    if ( nbnode == 1 ) {
-      if ( ID >= 0 ) e = mesh->Add0DElementWithID(node[0], ID);
-      else      e = mesh->Add0DElement      (node[0] );
-    }
-    break;
-  case SMDSAbs_Edge:
-    if ( nbnode == 2 ) {
-      if ( ID >= 0 ) e = mesh->AddEdgeWithID(node[0], node[1], ID);
-      else      e = mesh->AddEdge      (node[0], node[1] );
-    }
-    else if ( nbnode == 3 ) {
-      if ( ID >= 0 ) e = mesh->AddEdgeWithID(node[0], node[1], node[2], ID);
-      else      e = mesh->AddEdge      (node[0], node[1], node[2] );
-    }
-    break;
   case SMDSAbs_Face:
     if ( !isPoly ) {
       if      (nbnode == 3) {
-        if ( ID >= 0 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], ID);
-        else      e = mesh->AddFace      (node[0], node[1], node[2] );
+        if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], ID);
+        else           e = mesh->AddFace      (node[0], node[1], node[2] );
       }
       else if (nbnode == 4) {
-        if ( ID >= 0 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], ID);
-        else      e = mesh->AddFace      (node[0], node[1], node[2], node[3] );
+        if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], ID);
+        else           e = mesh->AddFace      (node[0], node[1], node[2], node[3] );
       }
       else if (nbnode == 6) {
-        if ( ID >= 0 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
-                                          node[4], node[5], ID);
-        else      e = mesh->AddFace      (node[0], node[1], node[2], node[3],
-                                          node[4], node[5] );
+        if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
+                                               node[4], node[5], ID);
+        else           e = mesh->AddFace      (node[0], node[1], node[2], node[3],
+                                               node[4], node[5] );
       }
       else if (nbnode == 8) {
-        if ( ID >= 0 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
-                                          node[4], node[5], node[6], node[7], ID);
-        else      e = mesh->AddFace      (node[0], node[1], node[2], node[3],
-                                          node[4], node[5], node[6], node[7] );
+        if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
+                                               node[4], node[5], node[6], node[7], ID);
+        else           e = mesh->AddFace      (node[0], node[1], node[2], node[3],
+                                               node[4], node[5], node[6], node[7] );
       }
     } else {
-      if ( ID >= 0 ) e = mesh->AddPolygonalFaceWithID(node, ID);
-      else      e = mesh->AddPolygonalFace      (node    );
+      if ( ID >= 1 ) e = mesh->AddPolygonalFaceWithID(node, ID);
+      else           e = mesh->AddPolygonalFace      (node    );
     }
     break;
+
   case SMDSAbs_Volume:
     if ( !isPoly ) {
       if      (nbnode == 4) {
-        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3] );
+        if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3] );
       }
       else if (nbnode == 5) {
-        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
-                                            node[4], ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
-                                            node[4] );
+        if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
+                                                 node[4], ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
+                                                 node[4] );
       }
       else if (nbnode == 6) {
-        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
-                                            node[4], node[5], ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
-                                            node[4], node[5] );
+        if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
+                                                 node[4], node[5], ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
+                                                 node[4], node[5] );
       }
       else if (nbnode == 8) {
-        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7], ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7] );
+        if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
+                                                 node[4], node[5], node[6], node[7], ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
+                                                 node[4], node[5], node[6], node[7] );
       }
       else if (nbnode == 10) {
-        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9], ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9] );
+        if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
+                                                 node[4], node[5], node[6], node[7],
+                                                 node[8], node[9], ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
+                                                 node[4], node[5], node[6], node[7],
+                                                 node[8], node[9] );
       }
       else if (nbnode == 13) {
-        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9], node[10],node[11],
-                                            node[12],ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9], node[10],node[11],
-                                            node[12] );
+        if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
+                                                 node[4], node[5], node[6], node[7],
+                                                 node[8], node[9], node[10],node[11],
+                                                 node[12],ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
+                                                 node[4], node[5], node[6], node[7],
+                                                 node[8], node[9], node[10],node[11],
+                                                 node[12] );
       }
       else if (nbnode == 15) {
-        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9], node[10],node[11],
-                                            node[12],node[13],node[14],ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9], node[10],node[11],
-                                            node[12],node[13],node[14] );
+        if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
+                                                 node[4], node[5], node[6], node[7],
+                                                 node[8], node[9], node[10],node[11],
+                                                 node[12],node[13],node[14],ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
+                                                 node[4], node[5], node[6], node[7],
+                                                 node[8], node[9], node[10],node[11],
+                                                 node[12],node[13],node[14] );
       }
       else if (nbnode == 20) {
-        if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9], node[10],node[11],
-                                            node[12],node[13],node[14],node[15],
-                                            node[16],node[17],node[18],node[19],ID);
-        else      e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
-                                            node[4], node[5], node[6], node[7],
-                                            node[8], node[9], node[10],node[11],
-                                            node[12],node[13],node[14],node[15],
-                                            node[16],node[17],node[18],node[19] );
+        if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
+                                                 node[4], node[5], node[6], node[7],
+                                                 node[8], node[9], node[10],node[11],
+                                                 node[12],node[13],node[14],node[15],
+                                                 node[16],node[17],node[18],node[19],ID);
+        else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
+                                                 node[4], node[5], node[6], node[7],
+                                                 node[8], node[9], node[10],node[11],
+                                                 node[12],node[13],node[14],node[15],
+                                                 node[16],node[17],node[18],node[19] );
       }
     }
+    break;
+
+  case SMDSAbs_Edge:
+    if ( nbnode == 2 ) {
+      if ( ID >= 1 ) e = mesh->AddEdgeWithID(node[0], node[1], ID);
+      else           e = mesh->AddEdge      (node[0], node[1] );
+    }
+    else if ( nbnode == 3 ) {
+      if ( ID >= 1 ) e = mesh->AddEdgeWithID(node[0], node[1], node[2], ID);
+      else           e = mesh->AddEdge      (node[0], node[1], node[2] );
+    }
+    break;
+
+  case SMDSAbs_0DElement:
+    if ( nbnode == 1 ) {
+      if ( ID >= 1 ) e = mesh->Add0DElementWithID(node[0], ID);
+      else           e = mesh->Add0DElement      (node[0] );
+    }
+    break;
+
+  case SMDSAbs_Node:
+    if ( ID >= 1 ) e = mesh->AddNodeWithID(node[0]->X(), node[0]->Y(), node[0]->Z(), ID);
+    else           e = mesh->AddNode      (node[0]->X(), node[0]->Y(), node[0]->Z());
+    break;
+
+  default:;
   }
   if ( e ) myLastCreatedElems.Append( e );
   return e;
@@ -352,45 +365,55 @@ int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem)
   if ( aMesh->ShapeToMesh().IsNull() )
     return 0;
 
-  if ( theElem->GetType() == SMDSAbs_Node )
-    {
-      int aShapeID = theElem->getshapeId();
-      if (aShapeID <= 0)
-        return 0;
-      else
-        return aShapeID;
-    }
-
-  TopoDS_Shape aShape; // the shape a node is on
-  SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
-  while ( nodeIt->more() ) {
-    const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
-    int aShapeID = node->getshapeId();
-    if (aShapeID > 0) {
-      SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID );
-      if ( sm ) {
-        if ( sm->Contains( theElem ))
-          return aShapeID;
-        if ( aShape.IsNull() )
-          aShape = aMesh->IndexToShape( aShapeID );
-      }
-      else {
-        //MESSAGE ( "::FindShape() No SubShape for aShapeID " << aShapeID );
+  int aShapeID = theElem->getshapeId();
+  if ( aShapeID < 1 )
+    return 0;
+
+  if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ))
+    if ( sm->Contains( theElem ))
+      return aShapeID;
+
+  if ( theElem->GetType() == SMDSAbs_Node ) {
+    MESSAGE( ":( Error: invalid myShapeId of node " << theElem->GetID() );
+  }
+  else {
+    MESSAGE( ":( Error: invalid myShapeId of element " << theElem->GetID() );
+  }
+
+  TopoDS_Shape aShape; // the shape a node of theElem is on
+  if ( theElem->GetType() != SMDSAbs_Node )
+  {
+    SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
+    while ( nodeIt->more() ) {
+      const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+      if ((aShapeID = node->getshapeId()) > 0) {
+        if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ) ) {
+          if ( sm->Contains( theElem ))
+            return aShapeID;
+          if ( aShape.IsNull() )
+            aShape = aMesh->IndexToShape( aShapeID );
+        }
       }
     }
   }
 
   // None of nodes is on a proper shape,
   // find the shape among ancestors of aShape on which a node is
-  if ( aShape.IsNull() ) {
-    //MESSAGE ("::FindShape() - NONE node is on shape")
-    return 0;
+  if ( !aShape.IsNull() ) {
+    TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape ));
+    for ( ; ancIt.More(); ancIt.Next() ) {
+      SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() );
+      if ( sm && sm->Contains( theElem ))
+        return aMesh->ShapeToIndex( ancIt.Value() );
+    }
   }
-  TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape ));
-  for ( ; ancIt.More(); ancIt.Next() ) {
-    SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() );
-    if ( sm && sm->Contains( theElem ))
-      return aMesh->ShapeToIndex( ancIt.Value() );
+  else
+  {
+    const map<int,SMESHDS_SubMesh*>& id2sm = GetMeshDS()->SubMeshes();
+    map<int,SMESHDS_SubMesh*>::const_iterator id_sm = id2sm.begin();
+    for ( ; id_sm != id2sm.end(); ++id_sm )
+      if ( id_sm->second->Contains( theElem ))
+        return id_sm->first;
   }
 
   //MESSAGE ("::FindShape() - SHAPE NOT FOUND")
@@ -1552,14 +1575,14 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
                                               const int                theMethodFlags)
 {
   // std-like iterator on coordinates of nodes of mesh element
-  typedef SMDS_StdIterator< TNodeXYZ, SMDS_ElemIteratorPtr > NXyzIterator;
+  typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > NXyzIterator;
   NXyzIterator xyzEnd;
 
   SMDS_VolumeTool    volTool;
   SMESH_MesherHelper helper( *GetMesh());
 
-  SMESHDS_SubMesh* subMesh = GetMeshDS()->MeshElements(1);
-  SMESHDS_SubMesh* fSubMesh = subMesh;
+  SMESHDS_SubMesh* subMesh = 0;//GetMeshDS()->MeshElements(1);
+  SMESHDS_SubMesh* fSubMesh = 0;//subMesh;
   
   SMESH_SequenceOfElemPtr newNodes, newElems;
 
@@ -1717,10 +1740,10 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
         }
         for ( int i = 0; i < triangles.size(); ++i )
         {
-          if ( !triangles.back() ) continue;
+          if ( !triangles[i] ) continue;
           if ( fSubMesh )
-            fSubMesh->AddElement( triangles.back());
-          newElems.Append( triangles.back() );
+            fSubMesh->AddElement( triangles[i]);
+          newElems.Append( triangles[i] );
         }
         ReplaceElemInGroups( face, triangles, GetMeshDS() );
         GetMeshDS()->RemoveFreeElement( face, fSubMesh, /*fromGroups=*/false );
@@ -2670,6 +2693,9 @@ void SMESH_MeshEditor::GetLinkedNodes( const SMDS_MeshNode* theNode,
   while ( elemIt->more() )
   {
     const SMDS_MeshElement* elem = elemIt->next();
+    if(elem->GetType() == SMDSAbs_0DElement)
+      continue;
+    
     SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
     if ( elem->GetType() == SMDSAbs_Volume )
     {
@@ -3874,6 +3900,8 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
     const SMDS_MeshElement* elem = itElem->first;
     vector<TNodeOfNodeListMapItr>& vecNewNodes = itElemNodes->second;
 
+    if(itElem->second.size()==0) continue;
+
     if ( elem->GetType() == SMDSAbs_Edge ) {
       // create a ceiling edge
       if (!elem->IsQuadratic()) {
@@ -3898,8 +3926,6 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
     if ( elem->GetType() != SMDSAbs_Face )
       continue;
 
-    if(itElem->second.size()==0) continue;
-
     bool hasFreeLinks = false;
 
     TIDSortedElemSet avoidSet;
@@ -5292,7 +5318,7 @@ void SMESH_MeshEditor::LinearAngleVariation(const int nbSteps,
  *  \param theCopy - if true, create translated copies of theElems
  *  \param theMakeGroups - if true and theCopy, create translated groups
  *  \param theTargetMesh - mesh to copy translated elements into
- *  \retval SMESH_MeshEditor::PGroupIDs - list of ids of created groups
+ *  \return SMESH_MeshEditor::PGroupIDs - list of ids of created groups
  */
 //================================================================================
 
@@ -6230,7 +6256,7 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
     const SMDS_MeshNode* closestNode = 0;
     list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
     for ( ; nIt != nodes.end(); ++nIt ) {
-      double sqDist = thePnt.SquareDistance( SMESH_MeshEditor::TNodeXYZ( *nIt ) );
+      double sqDist = thePnt.SquareDistance( SMESH_TNodeXYZ( *nIt ) );
       if ( minSqDist > sqDist ) {
         closestNode = *nIt;
         minSqDist = sqDist;
@@ -6285,7 +6311,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
   {
   public:
 
-    ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, double tolerance = NodeRadius );
+    ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr(), double tolerance = NodeRadius );
     void getElementsNearPoint( const gp_Pnt& point, TIDSortedElemSet& foundElems);
     void getElementsNearLine ( const gp_Ax1& line, TIDSortedElemSet& foundElems);
     ~ElementBndBoxTree();
@@ -6312,13 +6338,13 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
    */
   //================================================================================
 
-  ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, double tolerance)
+  ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, SMDS_ElemIteratorPtr theElemIt, double tolerance)
     :SMESH_Octree( new SMESH_Octree::Limit( MaxLevel, /*minSize=*/0. ))
   {
     int nbElems = mesh.GetMeshInfo().NbElements( elemType );
     _elements.reserve( nbElems );
 
-    SMDS_ElemIteratorPtr elemIt = mesh.elementsIterator( elemType );
+    SMDS_ElemIteratorPtr elemIt = theElemIt ? theElemIt : mesh.elementsIterator( elemType );
     while ( elemIt->more() )
       _elements.push_back( new ElementBox( elemIt->next(),tolerance  ));
 
@@ -6450,7 +6476,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
     _refCount = 1;
     SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
     while ( nIt->more() )
-      Add( SMESH_MeshEditor::TNodeXYZ( cast2Node( nIt->next() )));
+      Add( SMESH_TNodeXYZ( cast2Node( nIt->next() )));
     Enlarge( tolerance );
   }
 
@@ -6466,6 +6492,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint()
 struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
 {
   SMESHDS_Mesh*                _mesh;
+  SMDS_ElemIteratorPtr         _meshPartIt;
   ElementBndBoxTree*           _ebbTree;
   SMESH_NodeSearcherImpl*      _nodeSearcher;
   SMDSAbs_ElementType          _elementType;
@@ -6473,8 +6500,8 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
   bool                         _outerFacesFound;
   set<const SMDS_MeshElement*> _outerFaces; // empty means "no internal faces at all"
 
-  SMESH_ElementSearcherImpl( SMESHDS_Mesh& mesh )
-    : _mesh(&mesh),_ebbTree(0),_nodeSearcher(0), _tolerance(-1), _outerFacesFound(false) {}
+  SMESH_ElementSearcherImpl( SMESHDS_Mesh& mesh, SMDS_ElemIteratorPtr elemIt=SMDS_ElemIteratorPtr())
+    : _mesh(&mesh),_meshPartIt(elemIt),_ebbTree(0),_nodeSearcher(0),_tolerance(-1),_outerFacesFound(false) {}
   ~SMESH_ElementSearcherImpl()
   {
     if ( _ebbTree )      delete _ebbTree;      _ebbTree      = 0;
@@ -6556,7 +6583,7 @@ double SMESH_ElementSearcherImpl::getTolerance()
         SMDS_NodeIteratorPtr nodeIt = _mesh->nodesIterator();
         elemSize = 1;
         if ( meshInfo.NbNodes() > 2 )
-          elemSize = SMESH_MeshEditor::TNodeXYZ( nodeIt->next() ).Distance( nodeIt->next() );
+          elemSize = SMESH_TNodeXYZ( nodeIt->next() ).Distance( nodeIt->next() );
       }
       else
       {
@@ -6564,7 +6591,8 @@ double SMESH_ElementSearcherImpl::getTolerance()
             _mesh->elementsIterator( SMDSAbs_ElementType( complexType ));
         const SMDS_MeshElement* elem = elemIt->next();
         SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
-        SMESH_MeshEditor::TNodeXYZ n1( cast2Node( nodeIt->next() ));
+        SMESH_TNodeXYZ n1( cast2Node( nodeIt->next() ));
+        elemSize = 0;
         while ( nodeIt->more() )
         {
           double dist = n1.Distance( cast2Node( nodeIt->next() ));
@@ -6597,8 +6625,8 @@ bool SMESH_ElementSearcherImpl::getIntersParamOnLine(const gp_Lin&           lin
   int nbNodes = face->IsQuadratic() ? face->NbNodes()/2 : face->NbNodes();
   for ( int i = 0; i < nbNodes && nbInts < 2; ++i )
   {
-    GC_MakeSegment edge( SMESH_MeshEditor::TNodeXYZ( face->GetNode( i )),
-                         SMESH_MeshEditor::TNodeXYZ( face->GetNode( (i+1)%nbNodes) )); 
+    GC_MakeSegment edge( SMESH_TNodeXYZ( face->GetNode( i )),
+                         SMESH_TNodeXYZ( face->GetNode( (i+1)%nbNodes) )); 
     anExtCC.Init( lineCurve, edge);
     if ( anExtCC.NbExtrema() > 0 && anExtCC.LowerDistance() <= tol)
     {
@@ -6658,8 +6686,8 @@ void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerF
       seamLinks.insert( link );
 
       // link direction within the outerFace
-      gp_Vec n1n2( SMESH_MeshEditor::TNodeXYZ( link.node1()),
-                   SMESH_MeshEditor::TNodeXYZ( link.node2()));
+      gp_Vec n1n2( SMESH_TNodeXYZ( link.node1()),
+                   SMESH_TNodeXYZ( link.node2()));
       int i1 = outerFace->GetNodeIndex( link.node1() );
       int i2 = outerFace->GetNodeIndex( link.node2() );
       bool rev = ( abs(i2-i1) == 1 ? i1 > i2 : i2 > i1 );
@@ -6742,7 +6770,7 @@ FindElementsByPoint(const gp_Pnt&                      point,
     const SMDS_MeshNode* closeNode = _nodeSearcher->FindClosestTo( point );
     if ( !closeNode ) return foundElements.size();
 
-    if ( point.Distance( SMESH_MeshEditor::TNodeXYZ( closeNode )) > tolerance )
+    if ( point.Distance( SMESH_TNodeXYZ( closeNode )) > tolerance )
       return foundElements.size(); // to far from any node
 
     if ( type == SMDSAbs_Node )
@@ -6762,7 +6790,7 @@ FindElementsByPoint(const gp_Pnt&                      point,
     if ( !_ebbTree || _elementType != type )
     {
       if ( _ebbTree ) delete _ebbTree;
-      _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, tolerance );
+      _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt, tolerance );
     }
     TIDSortedElemSet suspectElems;
     _ebbTree->getElementsNearPoint( point, suspectElems );
@@ -6786,7 +6814,7 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
   if ( !_ebbTree || _elementType != SMDSAbs_Face )
   {
     if ( _ebbTree ) delete _ebbTree;
-    _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = SMDSAbs_Face );
+    _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = SMDSAbs_Face, _meshPartIt );
   }
   // Algo: analyse transition of a line starting at the point through mesh boundary;
   // try three lines parallel to axis of the coordinate system and perform rough
@@ -6814,7 +6842,7 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point)
       // get face plane
       gp_XYZ fNorm;
       if ( !SMESH_Algo::FaceNormal( *face, fNorm, /*normalized=*/false)) continue;
-      gp_Pln facePlane( SMESH_MeshEditor::TNodeXYZ( (*face)->GetNode(0)), fNorm );
+      gp_Pln facePlane( SMESH_TNodeXYZ( (*face)->GetNode(0)), fNorm );
 
       // perform intersection
       IntAna_IntConicQuad intersection( line, IntAna_Quadric( facePlane ));
@@ -7010,7 +7038,7 @@ void SMESH_ElementSearcherImpl::GetElementsNearLine( const gp_Ax1&
   if ( !_ebbTree || _elementType != type )
   {
     if ( _ebbTree ) delete _ebbTree;
-    _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type );
+    _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt );
   }
   TIDSortedElemSet suspectFaces; // elements possibly intersecting the line
   _ebbTree->getElementsNearLine( line, suspectFaces );
@@ -7028,6 +7056,17 @@ SMESH_ElementSearcher* SMESH_MeshEditor::GetElementSearcher()
   return new SMESH_ElementSearcherImpl( *GetMeshDS() );
 }
 
+//=======================================================================
+/*!
+ * \brief Return SMESH_ElementSearcher
+ */
+//=======================================================================
+
+SMESH_ElementSearcher* SMESH_MeshEditor::GetElementSearcher(SMDS_ElemIteratorPtr elemIt)
+{
+  return new SMESH_ElementSearcherImpl( *GetMeshDS(), elemIt );
+}
+
 //=======================================================================
 /*!
  * \brief Return true if the point is IN or ON of the element
@@ -7056,7 +7095,7 @@ bool SMESH_MeshEditor::isOut( const SMDS_MeshElement* element, const gp_Pnt& poi
   while ( nodeIt->more() )
     {
       const SMDS_MeshNode* node = cast2Node( nodeIt->next() );
-      xyz.push_back( TNodeXYZ(node) );
+      xyz.push_back( SMESH_TNodeXYZ(node) );
       nodeList.push_back(node);
     }
 
@@ -9156,7 +9195,7 @@ void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode*        theBetweenNode
 //=======================================================================
 /*!
  * \brief Convert elements contained in a submesh to quadratic
- * \retval int - nb of checked elements
+ * \return int - nb of checked elements
  */
 //=======================================================================
 
@@ -9176,8 +9215,6 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
     if( !elem || elem->IsQuadratic() ) continue;
 
     int id = elem->GetID();
-    //MESSAGE("elem " << id);
-    id = 0; // get a free number for new elements
     int nbNodes = elem->NbNodes();
     SMDSAbs_ElementType aType = elem->GetType();
 
@@ -9185,6 +9222,8 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
     if ( elem->GetEntityType() == SMDSEntity_Polyhedra )
       nbNodeInFaces = static_cast<const SMDS_VtkVolume* >( elem )->GetQuantities();
 
+    GetMeshDS()->RemoveFreeElement(elem, theSm, /*fromGroups=*/false);
+
     const SMDS_MeshElement* NewElem = 0;
 
     switch( aType )
@@ -9238,8 +9277,6 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh *   theSm,
     ReplaceElemInGroups( elem, NewElem, GetMeshDS());
     if( NewElem )
       theSm->AddElement( NewElem );
-
-    GetMeshDS()->RemoveFreeElement(elem, theSm, /*fromGroups=*/false);
   }
 //  if (!GetMeshDS()->isCompacted())
 //    GetMeshDS()->compactMesh();
@@ -9365,14 +9402,152 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
     aHelper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh
     aHelper.FixQuadraticElements();
   }
-  if (!GetMeshDS()->isCompacted())
-    GetMeshDS()->compactMesh();
+}
+
+//================================================================================
+/*!
+ * \brief Makes given elements quadratic
+ *  \param theForce3d - if true, the medium nodes will be placed in the middle of link
+ *  \param theElements - elements to make quadratic 
+ */
+//================================================================================
+
+void SMESH_MeshEditor::ConvertToQuadratic(const bool        theForce3d,
+                                          TIDSortedElemSet& theElements)
+{
+  if ( theElements.empty() ) return;
+
+  // we believe that all theElements are of the same type
+  SMDSAbs_ElementType elemType = (*theElements.begin())->GetType();
+  
+  // get all nodes shared by theElements
+  TIDSortedNodeSet allNodes;
+  TIDSortedElemSet::iterator eIt = theElements.begin();
+  for ( ; eIt != theElements.end(); ++eIt )
+    allNodes.insert( (*eIt)->begin_nodes(), (*eIt)->end_nodes() );
+
+  // complete theElements with elements of lower dim whose all nodes are in allNodes
+
+  TIDSortedElemSet quadAdjacentElems    [ SMDSAbs_NbElementTypes ]; // quadratic adjacent elements
+  TIDSortedElemSet checkedAdjacentElems [ SMDSAbs_NbElementTypes ];
+  TIDSortedNodeSet::iterator nIt = allNodes.begin();
+  for ( ; nIt != allNodes.end(); ++nIt )
+  {
+    const SMDS_MeshNode* n = *nIt;
+    SMDS_ElemIteratorPtr invIt = n->GetInverseElementIterator();
+    while ( invIt->more() )
+    {
+      const SMDS_MeshElement* e = invIt->next();
+      if ( e->IsQuadratic() )
+      {
+        quadAdjacentElems[ e->GetType() ].insert( e );
+        continue;
+      }
+      if ( e->GetType() >= elemType )
+      {
+        continue; // same type of more complex linear element
+      }
+
+      if ( !checkedAdjacentElems[ e->GetType() ].insert( e ).second )
+        continue; // e is already checked
+
+      // check nodes
+      bool allIn = true;
+      SMDS_ElemIteratorPtr nodeIt = e->nodesIterator();
+      while ( nodeIt->more() && allIn )
+        allIn = allNodes.count( cast2Node( nodeIt->next() ));
+      if ( allIn )
+        theElements.insert(e );
+    }
+  }
+
+  SMESH_MesherHelper helper(*myMesh);
+  helper.SetIsQuadratic( true );
+
+  // add links of quadratic adjacent elements to the helper
+
+  if ( !quadAdjacentElems[SMDSAbs_Edge].empty() )
+    for ( eIt  = quadAdjacentElems[SMDSAbs_Edge].begin();
+          eIt != quadAdjacentElems[SMDSAbs_Edge].end(); ++eIt )
+    {
+      helper.AddTLinks( static_cast< const SMDS_MeshEdge*> (*eIt) );
+    }
+  if ( !quadAdjacentElems[SMDSAbs_Face].empty() )
+    for ( eIt  = quadAdjacentElems[SMDSAbs_Face].begin();
+          eIt != quadAdjacentElems[SMDSAbs_Face].end(); ++eIt )
+    {
+      helper.AddTLinks( static_cast< const SMDS_MeshFace*> (*eIt) );
+    }
+  if ( !quadAdjacentElems[SMDSAbs_Volume].empty() )
+    for ( eIt  = quadAdjacentElems[SMDSAbs_Volume].begin();
+          eIt != quadAdjacentElems[SMDSAbs_Volume].end(); ++eIt )
+    {
+      helper.AddTLinks( static_cast< const SMDS_MeshVolume*> (*eIt) );
+    }
+
+  // make quadratic elements instead of linear ones
+
+  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() )
+      continue;
+
+    int id = elem->GetID();
+    SMDSAbs_ElementType type = elem->GetType();
+    vector<const SMDS_MeshNode *> nodes ( elem->begin_nodes(), elem->end_nodes());
+
+    if ( !smDS || !smDS->Contains( elem ))
+      smDS = meshDS->MeshElements( elem->getshapeId() );
+    meshDS->RemoveFreeElement(elem, smDS, /*fromGroups=*/false);
+
+    SMDS_MeshElement * newElem = 0;
+    switch( nodes.size() )
+    {
+    case 4: // cases for most multiple element types go first (for optimization)
+      if ( type == SMDSAbs_Volume )
+        newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+      else
+        newElem = helper.AddFace  (nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+      break;
+    case 8:
+      newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+                                 nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
+      break;
+    case 3:
+      newElem = helper.AddFace  (nodes[0], nodes[1], nodes[2], id, theForce3d);
+      break;
+    case 2:
+      newElem = helper.AddEdge(nodes[0], nodes[1], id, theForce3d);
+      break;
+    case 5:
+      newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+                                 nodes[4], id, theForce3d);
+      break;
+    case 6:
+      newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+                                 nodes[4], nodes[5], id, theForce3d);
+      break;
+    default:;
+    }
+    ReplaceElemInGroups( elem, newElem, meshDS);
+    if( newElem && smDS )
+      smDS->AddElement( newElem );
+  }
+
+  if ( !theForce3d  && !getenv("NO_FixQuadraticElements"))
+  { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
+    helper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh
+    helper.FixQuadraticElements();
+  }
 }
 
 //=======================================================================
 /*!
  * \brief Convert quadratic elements to linear ones and remove quadratic nodes
- * \retval int - nb of checked elements
+ * \return int - nb of checked elements
  */
 //=======================================================================
 
@@ -9382,7 +9557,6 @@ int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh *    theSm,
 {
   int nbElem = 0;
   SMESHDS_Mesh* meshDS = GetMeshDS();
-  const bool notFromGroups = false;
 
   while( theItr->more() )
   {
@@ -9390,44 +9564,28 @@ int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh *    theSm,
     nbElem++;
     if( elem && elem->IsQuadratic())
     {
-      int id = elem->GetID();
-      int nbNodes = elem->NbNodes();
-      vector<const SMDS_MeshNode *> nodes, mediumNodes;
-      nodes.reserve( nbNodes );
-      mediumNodes.reserve( nbNodes );
-
-      for(int i = 0; i < nbNodes; i++)
-      {
-        const SMDS_MeshNode* n = elem->GetNode(i);
-
-        if( elem->IsMediumNode( n ) )
-          mediumNodes.push_back( n );
-        else
-          nodes.push_back( n );
-      }
-      if( nodes.empty() ) continue;
+      int id                    = elem->GetID();
+      int nbCornerNodes         = elem->NbCornerNodes();
       SMDSAbs_ElementType aType = elem->GetType();
 
-      //remove old quadratic element
-      meshDS->RemoveFreeElement( elem, theSm, notFromGroups );
+      vector<const SMDS_MeshNode *> nodes( elem->begin_nodes(), elem->end_nodes() );
 
-      SMDS_MeshElement * NewElem = AddElement( nodes, aType, false, id );
-      ReplaceElemInGroups(elem, NewElem, meshDS);
-      if( theSm && NewElem )
-        theSm->AddElement( NewElem );
+      //remove a quadratic element
+      if ( !theSm || !theSm->Contains( elem ))
+        theSm = meshDS->MeshElements( elem->getshapeId() );
+      meshDS->RemoveFreeElement( elem, theSm, /*fromGroups=*/false );
 
       // remove medium nodes
-      vector<const SMDS_MeshNode*>::iterator nIt = mediumNodes.begin();
-      for ( ; nIt != mediumNodes.end(); ++nIt ) {
-        const SMDS_MeshNode* n = *nIt;
-        if ( n->NbInverseElements() == 0 ) {
-          if ( n->getshapeId() != theShapeID )
-            meshDS->RemoveFreeNode( n, meshDS->MeshElements
-                                    ( n->getshapeId() ));
-          else
-            meshDS->RemoveFreeNode( n, theSm );
-        }
-      }
+      for ( unsigned i = nbCornerNodes; i < nodes.size(); ++i )
+        if ( nodes[i]->NbInverseElements() == 0 )
+          meshDS->RemoveFreeNode( nodes[i], theSm );
+
+      // add a linear element
+      nodes.resize( nbCornerNodes );
+      SMDS_MeshElement * newElem = AddElement( nodes, aType, false, id );
+      ReplaceElemInGroups(elem, newElem, meshDS);
+      if( theSm && newElem )
+        theSm->AddElement( newElem );
     }
   }
   return nbElem;
@@ -9437,7 +9595,8 @@ int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh *    theSm,
 //function : ConvertFromQuadratic
 //purpose  :
 //=======================================================================
-bool  SMESH_MeshEditor::ConvertFromQuadratic()
+
+bool SMESH_MeshEditor::ConvertFromQuadratic()
 {
   int nbCheckedElems = 0;
   if ( myMesh->HasShapeToMesh() )
@@ -9464,6 +9623,102 @@ bool  SMESH_MeshEditor::ConvertFromQuadratic()
   return true;
 }
 
+namespace
+{
+  //================================================================================
+  /*!
+   * \brief Return true if all medium nodes of the element are in the node set
+   */
+  //================================================================================
+
+  bool allMediumNodesIn(const SMDS_MeshElement* elem, TIDSortedNodeSet& nodeSet )
+  {
+    for ( int i = elem->NbCornerNodes(); i < elem->NbNodes(); ++i )
+      if ( !nodeSet.count( elem->GetNode(i) ))
+        return false;
+    return true;
+  }
+}
+
+//================================================================================
+/*!
+ * \brief Makes given elements linear
+ */
+//================================================================================
+
+void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements)
+{
+  if ( theElements.empty() ) return;
+
+  // collect IDs of medium nodes of theElements; some of these nodes will be removed
+  set<int> mediumNodeIDs;
+  TIDSortedElemSet::iterator eIt = theElements.begin();
+  for ( ; eIt != theElements.end(); ++eIt )
+  {
+    const SMDS_MeshElement* e = *eIt;
+    for ( int i = e->NbCornerNodes(); i < e->NbNodes(); ++i )
+      mediumNodeIDs.insert( e->GetNode(i)->GetID() );
+  }
+
+  // replace given elements by linear ones
+  typedef SMDS_SetIterator<const SMDS_MeshElement*, TIDSortedElemSet::iterator> TSetIterator;
+  SMDS_ElemIteratorPtr elemIt( new TSetIterator( theElements.begin(), theElements.end() ));
+  removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 );
+
+  // we need to convert remaining elements whose all medium nodes are in mediumNodeIDs
+  // except those elements sharing medium nodes of quadratic element whose medium nodes
+  // are not all in mediumNodeIDs
+
+  // get remaining medium nodes
+  TIDSortedNodeSet mediumNodes;
+  set<int>::iterator nIdsIt = mediumNodeIDs.begin();
+  for ( ; nIdsIt != mediumNodeIDs.end(); ++nIdsIt )
+    if ( const SMDS_MeshNode* n = GetMeshDS()->FindNode( *nIdsIt ))
+      mediumNodes.insert( mediumNodes.end(), n );
+
+  // find more quadratic elements to convert
+  TIDSortedElemSet moreElemsToConvert;
+  TIDSortedNodeSet::iterator nIt = mediumNodes.begin();
+  for ( ; nIt != mediumNodes.end(); ++nIt )
+  {
+    SMDS_ElemIteratorPtr invIt = (*nIt)->GetInverseElementIterator();
+    while ( invIt->more() )
+    {
+      const SMDS_MeshElement* e = invIt->next();
+      if ( e->IsQuadratic() && allMediumNodesIn( e, mediumNodes ))
+      {
+        // find a more complex element including e and
+        // whose medium nodes are not in mediumNodes
+        bool complexFound = false;
+        for ( int type = e->GetType() + 1; type < SMDSAbs_0DElement; ++type )
+        {
+          SMDS_ElemIteratorPtr invIt2 =
+            (*nIt)->GetInverseElementIterator( SMDSAbs_ElementType( type ));
+          while ( invIt2->more() )
+          {
+            const SMDS_MeshElement* eComplex = invIt2->next();
+            if ( eComplex->IsQuadratic() && !allMediumNodesIn( eComplex, mediumNodes))
+            {
+              int nbCommonNodes = SMESH_Algo::GetCommonNodes( e, eComplex ).size();
+              if ( nbCommonNodes == e->NbNodes())
+              {
+                complexFound = true;
+                type = SMDSAbs_NbElementTypes; // to quit from the outer loop
+                break;
+              }
+            }
+          }
+        }
+        if ( !complexFound )
+          moreElemsToConvert.insert( e );
+      }
+    }
+  }
+  elemIt = SMDS_ElemIteratorPtr
+    (new TSetIterator( moreElemsToConvert.begin(), moreElemsToConvert.end() ));
+  removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 );
+}
+
 //=======================================================================
 //function : SewSideElements
 //purpose  :
@@ -10093,7 +10348,7 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet&    theSide1,
  * \param theSecondNode1 - a boundary node of set 1 linked with theFirstNode1
  * \param theSecondNode2 - a node of set 2 corresponding to theSecondNode1
  * \param nReplaceMap - output map of corresponding nodes
- * \retval bool  - is a success or not
+ * \return bool  - is a success or not
  */
 //================================================================================
 
@@ -10470,7 +10725,7 @@ namespace {
     gp_XYZ centerXYZ (0, 0, 0);
     SMDS_ElemIteratorPtr aNodeItr = theElem->nodesIterator();
     while (aNodeItr->more())
-      centerXYZ += SMESH_MeshEditor::TNodeXYZ(cast2Node( aNodeItr->next()));
+      centerXYZ += SMESH_TNodeXYZ(cast2Node( aNodeItr->next()));
 
     gp_Pnt aPnt = centerXYZ / theElem->NbNodes();
     theClassifier.Perform(aPnt, theTol);
@@ -10575,12 +10830,35 @@ bool SMESH_MeshEditor::DoubleNodesInRegion( const TIDSortedElemSet& theElems,
   return DoubleNodes( theElems, theNodesNot, anAffected );
 }
 
+/*!
+ *  \brief compute an oriented angle between two planes defined by four points.
+ *  The vector (p0,p1) defines the intersection of the 2 planes (p0,p1,g1) and (p0,p1,g2)
+ *  @param p0 base of the rotation axe
+ *  @param p1 extremity of the rotation axe
+ *  @param g1 belongs to the first plane
+ *  @param g2 belongs to the second plane
+ */
+double SMESH_MeshEditor::OrientedAngle(const gp_Pnt& p0, const gp_Pnt& p1, const gp_Pnt& g1, const gp_Pnt& g2)
+{
+//  MESSAGE("    p0: " << p0.X() << " " << p0.Y() << " " << p0.Z());
+//  MESSAGE("    p1: " << p1.X() << " " << p1.Y() << " " << p1.Z());
+//  MESSAGE("    g1: " << g1.X() << " " << g1.Y() << " " << g1.Z());
+//  MESSAGE("    g2: " << g2.X() << " " << g2.Y() << " " << g2.Z());
+  gp_Vec vref(p0, p1);
+  gp_Vec v1(p0, g1);
+  gp_Vec v2(p0, g2);
+  gp_Vec n1 = vref.Crossed(v1);
+  gp_Vec n2 = vref.Crossed(v2);
+  return n2.AngleWithRef(n1, vref);
+}
+
 /*!
  * \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.
  * @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
@@ -10589,23 +10867,29 @@ bool SMESH_MeshEditor::DoubleNodesInRegion( const TIDSortedElemSet& theElems,
 bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSortedElemSet>& theElems,
                                                      bool createJointElems)
 {
-  MESSAGE("------------------------------------------------------");
-  MESSAGE("SMESH_MeshEditor::CreateJointElementsOnGroupBoundaries");
-  MESSAGE("------------------------------------------------------");
+  MESSAGE("----------------------------------------------");
+  MESSAGE("SMESH_MeshEditor::doubleNodesOnGroupBoundaries");
+  MESSAGE("----------------------------------------------");
 
   SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS();
-  meshDS->BuildDownWardConnectivity(false);
+  meshDS->BuildDownWardConnectivity(true);
   CHRONO(50);
   SMDS_UnstructuredGrid *grid = meshDS->getGrid();
 
   // --- build the list of faces shared by 2 domains (group of elements), with their domain and volume indexes
+  //     build the list of cells with only a node or an edge on the border, with their domain and volume indexes
   //     build the list of nodes shared by 2 or more domains, with their domain indexes
 
-  std::map<DownIdType, std::map<int,int>, DownIdCompare> faceDomains; // 2x(id domain --> id volume)
-  std::map<int, std::map<int,int> > nodeDomains; //oldId ->  (domainId -> newId)
+  std::map<DownIdType, std::map<int,int>, DownIdCompare> faceDomains; // face --> (id domain --> id volume)
+  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();
 
   for (int idom = 0; idom < theElems.size(); idom++)
@@ -10640,43 +10924,217 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
                   if (!faceDomains[face].count(idom))
                     {
                       faceDomains[face][idom] = vtkId; // volume associated to face in this domain
+                      celldom[vtkId] = idom;
                     }
                 }
             }
         }
     }
 
-  MESSAGE("Number of shared faces " << faceDomains.size());
+  //MESSAGE("Number of shared faces " << faceDomains.size());
+  std::map<DownIdType, std::map<int, int>, DownIdCompare>::iterator itface;
 
-  // --- for each shared face, get the nodes
+  // --- explore the shared faces domain by domain,
+  //     explore the nodes of the face and see if they belong to a cell in the domain,
+  //     which has only a node or an edge on the border (not a shared face)
+
+  for (int idomain = 0; idomain < theElems.size(); idomain++)
+    {
+      const TIDSortedElemSet& domain = theElems[idomain];
+      itface = faceDomains.begin();
+      for (; itface != faceDomains.end(); ++itface)
+        {
+          std::map<int, int> domvol = itface->second;
+          if (!domvol.count(idomain))
+            continue;
+          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)
+            {
+              int oldId = *itn;
+              //MESSAGE("     node " << oldId);
+              std::set<int> cells;
+              cells.clear();
+              vtkCellLinks::Link l = grid->GetCellLinks()->GetLink(oldId);
+              for (int i=0; i<l.ncells; i++)
+                {
+                  int vtkId = l.cells[i];
+                  const SMDS_MeshElement* anElem = GetMeshDS()->FindElement(GetMeshDS()->fromVtkToSmds(vtkId));
+                  if (!domain.count(anElem))
+                    continue;
+                  int vtkType = grid->GetCellType(vtkId);
+                  int downId = grid->CellIdToDownId(vtkId);
+                  if (downId < 0)
+                    {
+                      MESSAGE("doubleNodesOnGroupBoundaries: internal algorithm problem");
+                      continue; // not OK at this stage of the algorithm:
+                                //no cells created after BuildDownWardConnectivity
+                    }
+                  DownIdType aCell(downId, vtkType);
+                  if (celldom.count(vtkId))
+                    continue;
+                  cellDomains[aCell][idomain] = vtkId;
+                  celldom[vtkId] = idomain;
+                }
+            }
+        }
+    }
+
+  // --- explore the shared faces domain by domain, to duplicate the nodes in a coherent way
+  //     for each shared face, get the nodes
   //     for each node, for each domain of the face, create a clone of the node
 
-  std::map<DownIdType, std::map<int,int>, DownIdCompare>::iterator itface = faceDomains.begin();
-  for( ; itface != faceDomains.end();++itface )
+  // --- edges at the intersection of 3 or 4 domains, with the order of domains to build
+  //     junction elements of type prism or hexa. the key is the pair of nodesId (lower first)
+  //     the value is the ordered domain ids. (more than 4 domains not taken into account)
+
+  std::map<std::vector<int>, std::vector<int> > edgesMultiDomains; // nodes of edge --> ordered domains
+  std::map<int, std::vector<int> > mutipleNodes; // nodes muti domains with domain order
+
+  for (int idomain = 0; idomain < theElems.size(); idomain++)
     {
-      DownIdType face = itface->first;
-      std::map<int,int> domvol = itface->second;
-      std::set<int> oldNodes;
-      oldNodes.clear();
-      grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
-      std::set<int>::iterator itn = oldNodes.begin();
-      for (;itn != oldNodes.end(); ++itn)
+      itface = faceDomains.begin();
+      for (; itface != faceDomains.end(); ++itface)
         {
-          int oldId = *itn;
-          if (!nodeDomains.count(oldId))
-            nodeDomains[oldId] = emptyMap; // create an empty entry for node
-          std::map<int,int>::iterator itdom = domvol.begin();
-          for(; itdom != domvol.end(); ++itdom)
+          std::map<int, int> domvol = itface->second;
+          if (!domvol.count(idomain))
+            continue;
+          DownIdType face = itface->first;
+          //MESSAGE(" --- face " << face.cellId);
+          std::set<int> oldNodes;
+          oldNodes.clear();
+          grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
+          bool isMultipleDetected = false;
+          std::set<int>::iterator itn = oldNodes.begin();
+          for (; itn != oldNodes.end(); ++itn)
             {
-              int idom = itdom->first;
-              if ( nodeDomains[oldId].empty() )
-                nodeDomains[oldId][idom] = oldId; // keep the old node in the first domain
-              else
+              int oldId = *itn;
+              //MESSAGE("     node " << oldId);
+              if (!nodeDomains.count(oldId))
+                nodeDomains[oldId] = emptyMap; // create an empty entry for node
+              if (nodeDomains[oldId].empty())
+                nodeDomains[oldId][idomain] = oldId; // keep the old node in the first domain
+              std::map<int, int>::iterator itdom = domvol.begin();
+              for (; itdom != domvol.end(); ++itdom)
+                {
+                  int idom = itdom->first;
+                  //MESSAGE("         domain " << idom);
+                  if (!nodeDomains[oldId].count(idom)) // --- node to clone
+                    {
+                      if (nodeDomains[oldId].size() >= 2) // a multiple node
+                        {
+                          vector<int> orderedDoms;
+                          //MESSAGE("multiple node " << oldId);
+                          isMultipleDetected =true;
+                          if (mutipleNodes.count(oldId))
+                            orderedDoms = mutipleNodes[oldId];
+                          else
+                            {
+                              map<int,int>::iterator it = nodeDomains[oldId].begin();
+                              for (; it != nodeDomains[oldId].end(); ++it)
+                                orderedDoms.push_back(it->first);
+                            }
+                          orderedDoms.push_back(idom); // TODO order ==> push_front or back
+                          //stringstream txt;
+                          //for (int i=0; i<orderedDoms.size(); i++)
+                          //  txt << orderedDoms[i] << " ";
+                          //MESSAGE("orderedDoms " << txt.str());
+                          mutipleNodes[oldId] = orderedDoms;
+                        }
+                      double *coords = grid->GetPoint(oldId);
+                      SMDS_MeshNode *newNode = meshDS->AddNode(coords[0], coords[1], coords[2]);
+                      int newId = newNode->getVtkId();
+                      nodeDomains[oldId][idom] = newId; // cloned node for other domains
+                      //MESSAGE("   newNode " << newId << " oldNode " << oldId << " size=" <<nodeDomains[oldId].size());
+                    }
+                  if (nodeDomains[oldId].size() >= 3)
+                    {
+                      //MESSAGE("confirm multiple node " << oldId);
+                      isMultipleDetected =true;
+                    }
+                }
+            }
+          if (isMultipleDetected) // check if an edge of the face is shared between 3 or more domains
+            {
+              //MESSAGE("multiple Nodes detected on a shared face");
+              int downId = itface->first.cellId;
+              unsigned char cellType = itface->first.cellType;
+              int nbEdges = grid->getDownArray(cellType)->getNumberOfDownCells(downId);
+              const int *downEdgeIds = grid->getDownArray(cellType)->getDownCells(downId);
+              const unsigned char* edgeType = grid->getDownArray(cellType)->getDownTypes(downId);
+              for (int ie =0; ie < nbEdges; ie++)
                 {
-                  double *coords = grid->GetPoint(oldId);
-                  SMDS_MeshNode *newNode = meshDS->AddNode(coords[0], coords[1], coords[2]);
-                  int newId = newNode->getVtkId();
-                  nodeDomains[oldId][idom] = newId; // cloned node for other domains
+                  int nodes[3];
+                  int nbNodes = grid->getDownArray(edgeType[ie])->getNodes(downEdgeIds[ie], nodes);
+                  if (mutipleNodes.count(nodes[0]) && mutipleNodes.count(nodes[nbNodes-1]))
+                    {
+                      vector<int> vn0 = mutipleNodes[nodes[0]];
+                      vector<int> vn1 = mutipleNodes[nodes[nbNodes - 1]];
+                      sort( vn0.begin(), vn0.end() );
+                      sort( vn1.begin(), vn1.end() );
+                      if (vn0 == vn1)
+                        {
+                          //MESSAGE(" detect edgesMultiDomains " << nodes[0] << " " << nodes[nbNodes - 1]);
+                          double *coords = grid->GetPoint(nodes[0]);
+                          gp_Pnt p0(coords[0], coords[1], coords[2]);
+                          coords = grid->GetPoint(nodes[nbNodes - 1]);
+                          gp_Pnt p1(coords[0], coords[1], coords[2]);
+                          gp_Pnt gref;
+                          int vtkVolIds[1000];  // an edge can belong to a lot of volumes
+                          map<int, SMDS_VtkVolume*> domvol; // domain --> a volume with the edge
+                          map<int, double> angleDom; // oriented angles between planes defined by edge and volume centers
+                          int nbvol = grid->GetParentVolumes(vtkVolIds, downEdgeIds[ie], edgeType[ie]);
+                          for (int id=0; id < vn0.size(); id++)
+                            {
+                              int idom = vn0[id];
+                              for (int ivol=0; ivol<nbvol; ivol++)
+                                {
+                                  int smdsId = meshDS->fromVtkToSmds(vtkVolIds[ivol]);
+                                  SMDS_MeshElement* elem = (SMDS_MeshElement*)meshDS->FindElement(smdsId);
+                                  if (theElems[idom].count(elem))
+                                    {
+                                      SMDS_VtkVolume* svol = dynamic_cast<SMDS_VtkVolume*>(elem);
+                                      domvol[idom] = svol;
+                                      //MESSAGE("  domain " << idom << " volume " << elem->GetID());
+                                      double values[3];
+                                      vtkIdType npts = 0;
+                                      vtkIdType* pts = 0;
+                                      grid->GetCellPoints(vtkVolIds[ivol], npts, pts);
+                                      SMDS_VtkVolume::gravityCenter(grid, pts, npts, values);
+                                      if (id ==0)
+                                        {
+                                          gref.SetXYZ(gp_XYZ(values[0], values[1], values[2]));
+                                          angleDom[idom] = 0;
+                                        }
+                                      else
+                                        {
+                                          gp_Pnt g(values[0], values[1], values[2]);
+                                          angleDom[idom] = OrientedAngle(p0, p1, gref, g); // -pi<angle<+pi
+                                          //MESSAGE("  angle=" << angleDom[idom]);
+                                        }
+                                      break;
+                                    }
+                                }
+                            }
+                          map<double, int> sortedDom; // sort domains by angle
+                          for (map<int, double>::iterator ia = angleDom.begin(); ia != angleDom.end(); ++ia)
+                            sortedDom[ia->second] = ia->first;
+                          vector<int> vnodes;
+                          vector<int> vdom;
+                          for (map<double, int>::iterator ib = sortedDom.begin(); ib != sortedDom.end(); ++ib)
+                            {
+                              vdom.push_back(ib->second);
+                              //MESSAGE("  ordered domain " << ib->second << "  angle " << ib->first);
+                            }
+                          for (int ino = 0; ino < nbNodes; ino++)
+                            vnodes.push_back(nodes[ino]);
+                          edgesMultiDomains[vnodes] = vdom; // nodes vector --> ordered domains
+                        }
+                    }
                 }
             }
         }
@@ -10686,42 +11144,128 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   //     get node id's of the face (id SMDS = id VTK)
   //     create flat element with old and new nodes if requested
 
+  // --- new quad nodes on flat quad elements: oldId --> ((domain1 X domain2) --> newId)
+  //     (domain1 X domain2) = domain1 + MAXINT*domain2
+
+  std::map<int, std::map<long,int> > nodeQuadDomains;
+  std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
+
   if (createJointElems)
     {
       itface = faceDomains.begin();
-      for( ; itface != faceDomains.end();++itface )
+      for (; itface != faceDomains.end(); ++itface)
         {
           DownIdType face = itface->first;
           std::set<int> oldNodes;
           std::set<int>::iterator itn;
           oldNodes.clear();
           grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
-          std::map<int,int> localClonedNodeIds;
 
-          std::map<int,int> domvol = itface->second;
-          std::map<int,int>::iterator itdom = domvol.begin();
+          std::map<int, int> domvol = itface->second;
+          std::map<int, int>::iterator itdom = domvol.begin();
           int dom1 = itdom->first;
           int vtkVolId = itdom->second;
           itdom++;
           int dom2 = itdom->first;
+          SMDS_MeshVolume *vol = grid->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains,
+                                                             nodeQuadDomains);
+          stringstream grpname;
+          grpname << "j_";
+          if (dom1 < dom2)
+            grpname << dom1 << "_" << dom2;
+          else
+            grpname << dom2 << "_" << dom1;
+          int idg;
+          string namegrp = grpname.str();
+          if (!mapOfJunctionGroups.count(namegrp))
+            mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg);
+          SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
+          if (sgrp)
+            sgrp->Add(vol->GetID());
+        }
+    }
+
+  // --- create volumes on multiple domain intersection if requested
+  //     iterate on edgesMultiDomains
 
-          localClonedNodeIds.clear();
-          for (itn = oldNodes.begin(); itn != oldNodes.end(); ++itn)
+  if (createJointElems)
+    {
+      std::map<std::vector<int>, std::vector<int> >::iterator ite = edgesMultiDomains.begin();
+      for (; ite != edgesMultiDomains.end(); ++ite)
+        {
+          vector<int> nodes = ite->first;
+          vector<int> orderDom = ite->second;
+          vector<vtkIdType> orderedNodes;
+          if (nodes.size() == 2)
             {
-              int oldId = *itn;
-              int refid = oldId;
-              if (nodeDomains[oldId].count(dom1))
-                refid = nodeDomains[oldId][dom1];
-              else
-                MESSAGE("--- problem domain node " << dom1 << " " << oldId);
-              int newid = oldId;
-              if (nodeDomains[oldId].count(dom2))
-                newid = nodeDomains[oldId][dom2];
-              else
-                MESSAGE("--- problem domain node " << dom2 << " " << oldId);
-              localClonedNodeIds[oldId] = newid;
+              //MESSAGE(" use edgesMultiDomains " << nodes[0] << " " << nodes[1]);
+              for (int ino=0; ino < nodes.size(); ino++)
+                if (orderDom.size() == 3)
+                  for (int idom = 0; idom <orderDom.size(); idom++)
+                    orderedNodes.push_back( nodeDomains[nodes[ino]][orderDom[idom]] );
+                else
+                  for (int idom = orderDom.size()-1; idom >=0; idom--)
+                    orderedNodes.push_back( nodeDomains[nodes[ino]][orderDom[idom]] );
+              SMDS_MeshVolume* vol = this->GetMeshDS()->AddVolumeFromVtkIds(orderedNodes);
+
+              stringstream grpname;
+              grpname << "mj_";
+              grpname << 0 << "_" << 0;
+              int idg;
+              string namegrp = grpname.str();
+              if (!mapOfJunctionGroups.count(namegrp))
+                mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg);
+              SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
+              if (sgrp)
+                sgrp->Add(vol->GetID());
+            }
+          else
+            {
+              //MESSAGE("Quadratic multiple joints not implemented");
+              // TODO quadratic nodes
+            }
+        }
+    }
+
+  // --- list the explicit faces and edges of the mesh that need to be modified,
+  //     i.e. faces and edges built with one or more duplicated nodes.
+  //     associate these faces or edges to their corresponding domain.
+  //     only the first domain found is kept when a face or edge is shared
+
+  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();
+
+  for (int idomain = 0; idomain < theElems.size(); idomain++)
+    {
+      std::map<int, std::map<int, int> >::const_iterator itnod = nodeDomains.begin();
+      for (; itnod != nodeDomains.end(); ++itnod)
+        {
+          int oldId = itnod->first;
+          //MESSAGE("     node " << oldId);
+          vtkCellLinks::Link l = grid->GetCellLinks()->GetLink(oldId);
+          for (int i = 0; i < l.ncells; i++)
+            {
+              int vtkId = l.cells[i];
+              int vtkType = grid->GetCellType(vtkId);
+              int downId = grid->CellIdToDownId(vtkId);
+              if (downId < 0)
+                continue; // new cells: not to be modified
+              DownIdType aCell(downId, vtkType);
+              int volParents[1000];
+              int nbvol = grid->GetParentVolumes(volParents, vtkId);
+              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);
+                    }
             }
-          meshDS->extrudeVolumeFromFace(vtkVolId, localClonedNodeIds);
         }
     }
 
@@ -10729,40 +11273,212 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector<TIDSorted
   //     get node id's of the face
   //     replace old nodes by new nodes in volumes, and update inverse connectivity
 
-  itface = faceDomains.begin();
-  for( ; itface != faceDomains.end();++itface )
+  std::map<DownIdType, std::map<int,int>, DownIdCompare>* maps[3] = {&faceDomains, &cellDomains, &faceOrEdgeDom};
+  for (int m=0; m<3; m++)
     {
-      DownIdType face = itface->first;
-      std::set<int> oldNodes;
-      std::set<int>::iterator itn;
-      oldNodes.clear();
-      grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
-      std::map<int,int> localClonedNodeIds;
-
-      std::map<int,int> domvol = itface->second;
-      std::map<int,int>::iterator itdom = domvol.begin();
-      for(; itdom != domvol.end(); ++itdom)
+      std::map<DownIdType, std::map<int,int>, DownIdCompare>* amap = maps[m];
+      itface = (*amap).begin();
+      for (; itface != (*amap).end(); ++itface)
         {
-          int idom = itdom->first;
-          int vtkVolId = itdom->second;
-          localClonedNodeIds.clear();
-          for (itn = oldNodes.begin(); itn != oldNodes.end(); ++itn)
+          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;
+
+          std::map<int, int> domvol = itface->second;
+          std::map<int, int>::iterator itdom = domvol.begin();
+          for (; itdom != domvol.end(); ++itdom)
             {
-              int oldId = *itn;
-              if (nodeDomains[oldId].count(idom))
-                localClonedNodeIds[oldId] = nodeDomains[oldId][idom];
+              int idom = itdom->first;
+              int vtkVolId = itdom->second;
+              //MESSAGE("modify nodes of cell " << this->GetMeshDS()->fromVtkToSmds(vtkVolId) << " domain " << idom);
+              localClonedNodeIds.clear();
+              for (itn = oldNodes.begin(); itn != oldNodes.end(); ++itn)
+                {
+                  int oldId = *itn;
+                  if (nodeDomains[oldId].count(idom))
+                    {
+                      localClonedNodeIds[oldId] = nodeDomains[oldId][idom];
+                      //MESSAGE("     node " << oldId << " --> " << localClonedNodeIds[oldId]);
+                    }
+                }
+              meshDS->ModifyCellNodes(vtkVolId, localClonedNodeIds);
             }
-          meshDS->ModifyCellNodes(vtkVolId, localClonedNodeIds);
         }
     }
+
+  meshDS->CleanDownWardConnectivity(); // Mesh has been modified, downward connectivity is no more usable, free memory
   grid->BuildLinks();
 
-  // TODO replace also old nodes by new nodes in faces and edges
   CHRONOSTOP(50);
   counters::stats();
   return true;
 }
 
+/*!
+ * \brief Double nodes on some external faces and create flat elements.
+ * Flat elements are mainly used by some types of mechanic calculations.
+ *
+ * Each group of the list must be constituted of faces.
+ * Triangles are transformed in prisms, and quadrangles in hexahedrons.
+ * @param theElems - list of groups of faces, where a group of faces is a set of
+ * SMDS_MeshElements sorted by Id.
+ * @return TRUE if operation has been completed successfully, FALSE otherwise
+ */
+bool SMESH_MeshEditor::CreateFlatElementsOnFacesGroups(const std::vector<TIDSortedElemSet>& theElems)
+{
+  MESSAGE("-------------------------------------------------");
+  MESSAGE("SMESH_MeshEditor::CreateFlatElementsOnFacesGroups");
+  MESSAGE("-------------------------------------------------");
+
+  SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS();
+
+  // --- For each group of faces
+  //     duplicate the nodes, create a flat element based on the face
+  //     replace the nodes of the faces by their clones
+
+  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 (int idom = 0; idom < theElems.size(); idom++)
+    {
+      const TIDSortedElemSet& domain = theElems[idom];
+      TIDSortedElemSet::const_iterator elemItr = domain.begin();
+      for (; elemItr != domain.end(); ++elemItr)
+        {
+          SMDS_MeshElement* anElem = (SMDS_MeshElement*) *elemItr;
+          SMDS_MeshFace* aFace = dynamic_cast<SMDS_MeshFace*> (anElem);
+          if (!aFace)
+            continue;
+          // MESSAGE("aFace=" << aFace->GetID());
+          bool isQuad = aFace->IsQuadratic();
+          vector<const SMDS_MeshNode*> ln0, ln1, ln2, ln3, ln4;
+
+          // --- clone the nodes, create intermediate nodes for non medium nodes of a quad face
+
+          SMDS_ElemIteratorPtr nodeIt = aFace->nodesIterator();
+          while (nodeIt->more())
+            {
+              const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*> (nodeIt->next());
+              bool isMedium = isQuad && (aFace->IsMediumNode(node));
+              if (isMedium)
+                ln2.push_back(node);
+              else
+                ln0.push_back(node);
+
+              const SMDS_MeshNode* clone = 0;
+              if (!clonedNodes.count(node))
+                {
+                  clone = meshDS->AddNode(node->X(), node->Y(), node->Z());
+                  clonedNodes[node] = clone;
+                }
+              else
+                clone = clonedNodes[node];
+
+              if (isMedium)
+                ln3.push_back(clone);
+              else
+                ln1.push_back(clone);
+
+              const SMDS_MeshNode* inter = 0;
+              if (isQuad && (!isMedium))
+                {
+                  if (!intermediateNodes.count(node))
+                    {
+                      inter = meshDS->AddNode(node->X(), node->Y(), node->Z());
+                      intermediateNodes[node] = inter;
+                    }
+                  else
+                    inter = intermediateNodes[node];
+                  ln4.push_back(inter);
+                }
+            }
+
+          // --- extrude the face
+
+          vector<const SMDS_MeshNode*> ln;
+          SMDS_MeshVolume* vol = 0;
+          vtkIdType aType = aFace->GetVtkType();
+          switch (aType)
+          {
+            case VTK_TRIANGLE:
+              vol = meshDS->AddVolume(ln0[2], ln0[1], ln0[0], ln1[2], ln1[1], ln1[0]);
+              // MESSAGE("vol prism " << vol->GetID());
+              ln.push_back(ln1[0]);
+              ln.push_back(ln1[1]);
+              ln.push_back(ln1[2]);
+              break;
+            case VTK_QUAD:
+              vol = meshDS->AddVolume(ln0[3], ln0[2], ln0[1], ln0[0], ln1[3], ln1[2], ln1[1], ln1[0]);
+              // MESSAGE("vol hexa " << vol->GetID());
+              ln.push_back(ln1[0]);
+              ln.push_back(ln1[1]);
+              ln.push_back(ln1[2]);
+              ln.push_back(ln1[3]);
+              break;
+            case VTK_QUADRATIC_TRIANGLE:
+              vol = meshDS->AddVolume(ln1[0], ln1[1], ln1[2], ln0[0], ln0[1], ln0[2], ln3[0], ln3[1], ln3[2],
+                                      ln2[0], ln2[1], ln2[2], ln4[0], ln4[1], ln4[2]);
+              // MESSAGE("vol quad prism " << vol->GetID());
+              ln.push_back(ln1[0]);
+              ln.push_back(ln1[1]);
+              ln.push_back(ln1[2]);
+              ln.push_back(ln3[0]);
+              ln.push_back(ln3[1]);
+              ln.push_back(ln3[2]);
+              break;
+            case VTK_QUADRATIC_QUAD:
+//              vol = meshDS->AddVolume(ln0[0], ln0[1], ln0[2], ln0[3], ln1[0], ln1[1], ln1[2], ln1[3],
+//                                      ln2[0], ln2[1], ln2[2], ln2[3], ln3[0], ln3[1], ln3[2], ln3[3],
+//                                      ln4[0], ln4[1], ln4[2], ln4[3]);
+              vol = meshDS->AddVolume(ln1[0], ln1[1], ln1[2], ln1[3], ln0[0], ln0[1], ln0[2], ln0[3],
+                                      ln3[0], ln3[1], ln3[2], ln3[3], ln2[0], ln2[1], ln2[2], ln2[3],
+                                      ln4[0], ln4[1], ln4[2], ln4[3]);
+              // MESSAGE("vol quad hexa " << vol->GetID());
+              ln.push_back(ln1[0]);
+              ln.push_back(ln1[1]);
+              ln.push_back(ln1[2]);
+              ln.push_back(ln1[3]);
+              ln.push_back(ln3[0]);
+              ln.push_back(ln3[1]);
+              ln.push_back(ln3[2]);
+              ln.push_back(ln3[3]);
+              break;
+            case VTK_POLYGON:
+              break;
+            default:
+              break;
+          }
+
+          if (vol)
+            {
+              stringstream grpname;
+              grpname << "jf_";
+              grpname << idom;
+              int idg;
+              string namegrp = grpname.str();
+              if (!mapOfJunctionGroups.count(namegrp))
+                mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg);
+              SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
+              if (sgrp)
+                sgrp->Add(vol->GetID());
+            }
+
+          // --- modify the face
+
+          aFace->ChangeNodes(&ln[0], ln.size());
+        }
+    }
+  return true;
+}
+
 //================================================================================
 /*!
  * \brief Generates skin mesh (containing 2D cells) from 3D mesh
@@ -10833,15 +11549,23 @@ namespace
  *  \param toCopyElements - if true, the checked elements will be copied into the targetMesh
  *  \param toCopyExistingBondary - if true, not only new but also pre-existing
  *                                boundary elements will be copied into the targetMesh
+ *  \param toAddExistingBondary - if true, not only new but also pre-existing
+ *                                boundary elements will be added into the new group
+ *  \param aroundElements - if true, elements will be created on boundary of given
+ *                          elements else, on boundary of the whole mesh. This
+ *                          option works for 2D elements only.
+ * \return nb of added boundary elements
  */
 //================================================================================
 
-void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
-                                        Bnd_Dimension           dimension,
-                                        SMESH_Group*            group/*=0*/,
-                                        SMESH_Mesh*             targetMesh/*=0*/,
-                                        bool                    toCopyElements/*=false*/,
-                                        bool                    toCopyExistingBondary/*=false*/)
+int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
+                                       Bnd_Dimension           dimension,
+                                       SMESH_Group*            group/*=0*/,
+                                       SMESH_Mesh*             targetMesh/*=0*/,
+                                       bool                    toCopyElements/*=false*/,
+                                       bool                    toCopyExistingBondary/*=false*/,
+                                       bool                    toAddExistingBondary/*= false*/,
+                                       bool                    aroundElements/*= false*/)
 {
   SMDSAbs_ElementType missType = (dimension == BND_2DFROM3D) ? SMDSAbs_Face : SMDSAbs_Edge;
   SMDSAbs_ElementType elemType = (dimension == BND_1DFROM2D) ? SMDSAbs_Face : SMDSAbs_Volume;
@@ -10849,14 +11573,24 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
   if ( !elements.empty() && (*elements.begin())->GetType() != elemType )
     throw SALOME_Exception(LOCALIZED("wrong element type"));
 
+  if ( aroundElements && elemType == SMDSAbs_Volume )
+    throw SALOME_Exception(LOCALIZED("wrong element type for aroundElements==true"));
+
   if ( !targetMesh )
     toCopyElements = toCopyExistingBondary = false;
 
   SMESH_MeshEditor tgtEditor( targetMesh ? targetMesh : myMesh );
   SMESHDS_Mesh* aMesh = GetMeshDS(), *tgtMeshDS = tgtEditor.GetMeshDS();
+  int nbAddedBnd = 0;
+
+  // editor adding present bnd elements and optionally holding elements to add to the group
+  SMESH_MeshEditor* presentEditor;
+  SMESH_MeshEditor tgtEditor2( tgtEditor.GetMesh() );
+  presentEditor = toAddExistingBondary ? &tgtEditor : &tgtEditor2;
 
   SMDS_VolumeTool vTool;
-  TIDSortedElemSet emptySet, avoidSet;
+  TIDSortedElemSet avoidSet;
+  const TIDSortedElemSet emptySet, *elemSet = aroundElements ? &elements : &emptySet;
   int inode;
 
   typedef vector<const SMDS_MeshNode*> TConnectivity;
@@ -10872,7 +11606,9 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
     const SMDS_MeshElement* elem = eIt->next();
     const int iQuad = elem->IsQuadratic();
 
+    // ------------------------------------------------------------------------------------
     // 1. For an elem, get present bnd elements and connectivities of missing bnd elements
+    // ------------------------------------------------------------------------------------
     vector<const SMDS_MeshElement*> presentBndElems;
     vector<TConnectivity>           missingBndElems;
     TConnectivity nodes;
@@ -10924,7 +11660,7 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
       {
         nodes[0] = elem->GetNode(i);
         nodes[1] = elem->GetNode((i+1)%nbNodes);
-        if ( FindFaceInSet( nodes[0], nodes[1], emptySet, avoidSet))
+        if ( FindFaceInSet( nodes[0], nodes[1], *elemSet, avoidSet))
           continue; // not free link
 
         //if ( iQuad )
@@ -10937,7 +11673,9 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
       }
     }
 
+    // ---------------------------------
     // 2. Add missing boundary elements
+    // ---------------------------------
     if ( targetMesh != myMesh )
       // instead of making a map of nodes in this mesh and targetMesh,
       // we create nodes with same IDs. We can renumber them later, if needed
@@ -10947,16 +11685,28 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
         TConnectivity  nodes( srcNodes.size() );
         for ( inode = 0; inode < nodes.size(); ++inode )
           nodes[inode] = getNodeWithSameID( tgtMeshDS, srcNodes[inode] );
+        if ( aroundElements && tgtEditor.GetMeshDS()->FindElement( nodes,
+                                                                   missType,
+                                                                   /*noMedium=*/true))
+          continue;
         tgtEditor.AddElement(nodes, missType, elem->IsPoly() && nodes.size()/(iQuad+1)>4);
+        ++nbAddedBnd;
       }
     else
       for ( int i = 0; i < missingBndElems.size(); ++i )
       {
-        TConnectivity&  nodes = missingBndElems[i];
+        TConnectivity& nodes = missingBndElems[i];
+        if ( aroundElements && tgtEditor.GetMeshDS()->FindElement( nodes,
+                                                                   missType,
+                                                                   /*noMedium=*/true))
+          continue;
         tgtEditor.AddElement(nodes, missType, elem->IsPoly() && nodes.size()/(iQuad+1)>4);
+        ++nbAddedBnd;
       }
 
+    // ----------------------------------
     // 3. Copy present boundary elements
+    // ----------------------------------
     if ( toCopyExistingBondary )
       for ( int i = 0 ; i < presentBndElems.size(); ++i )
       {
@@ -10964,13 +11714,19 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
         TConnectivity nodes( e->NbNodes() );
         for ( inode = 0; inode < nodes.size(); ++inode )
           nodes[inode] = getNodeWithSameID( tgtMeshDS, e->GetNode(inode) );
-        tgtEditor.AddElement(nodes, missType, e->IsPoly());
-        // leave only missing elements in tgtEditor.myLastCreatedElems
-        tgtEditor.myLastCreatedElems.Remove( tgtEditor.myLastCreatedElems.Size() );
+        presentEditor->AddElement(nodes, missType, e->IsPoly());
+      }
+    else // store present elements to add them to a group
+      for ( int i = 0 ; i < presentBndElems.size(); ++i )
+      {
+        presentEditor->myLastCreatedElems.Append(presentBndElems[i]);
       }
+      
   } // loop on given elements
 
-  // 4. Fill group with missing boundary elements
+  // ---------------------------------------------
+  // 4. Fill group with boundary elements
+  // ---------------------------------------------
   if ( group )
   {
     if ( SMESHDS_Group* g = dynamic_cast<SMESHDS_Group*>( group->GetGroupDS() ))
@@ -10978,9 +11734,12 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
         g->SMDSGroup().Add( tgtEditor.myLastCreatedElems( i+1 ));
   }
   tgtEditor.myLastCreatedElems.Clear();
+  tgtEditor2.myLastCreatedElems.Clear();
 
+  // -----------------------
   // 5. Copy given elements
-  if ( toCopyElements )
+  // -----------------------
+  if ( toCopyElements && targetMesh != myMesh )
   {
     if (elements.empty())
       eIt = aMesh->elementsIterator(elemType);
@@ -10997,5 +11756,5 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
       tgtEditor.myLastCreatedElems.Clear();
     }
   }
-  return;
+  return nbAddedBnd;
 }