diff --git a/CMakeLists.txt b/CMakeLists.txt
index ad268dcd0..9409ace73 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -27,7 +27,7 @@ INCLUDE(CMakeDependentOption)
 STRING(TOUPPER ${PROJECT_NAME} PROJECT_NAME_UC)
 
 SET(${PROJECT_NAME_UC}_MAJOR_VERSION 9)
-SET(${PROJECT_NAME_UC}_MINOR_VERSION 4)
+SET(${PROJECT_NAME_UC}_MINOR_VERSION 5)
 SET(${PROJECT_NAME_UC}_PATCH_VERSION 0)
 SET(${PROJECT_NAME_UC}_VERSION
   ${${PROJECT_NAME_UC}_MAJOR_VERSION}.${${PROJECT_NAME_UC}_MINOR_VERSION}.${${PROJECT_NAME_UC}_PATCH_VERSION})
diff --git a/src/SMDS/chrono.cxx b/src/SMDS/chrono.cxx
index ae198d95f..c4ad47ec4 100644
--- a/src/SMDS/chrono.cxx
+++ b/src/SMDS/chrono.cxx
@@ -21,6 +21,7 @@
 #include "utilities.h"
 
 using namespace std;
+using namespace SMDS;
 
 cntStruct* counters::_ctrs = 0;
 int counters::_nbChrono = 0;
diff --git a/src/SMDS/chrono.hxx b/src/SMDS/chrono.hxx
index 3fdae2eb2..032ac935c 100644
--- a/src/SMDS/chrono.hxx
+++ b/src/SMDS/chrono.hxx
@@ -27,6 +27,8 @@
 #include <iostream>
 #include <ctime>
 
+namespace SMDS
+{
 typedef struct acnt
 {
   char*  _ctrNames;
@@ -57,9 +59,11 @@ protected:
   int _ctr;
   clock_t _start, _end;
 };
-
+}
 #ifdef CHRONODEF
-#define CHRONO(i) counters::_ctrs[i]._ctrNames = (char *)__FILE__; \
+#define CHRONO(i)
+  using namespace SMDS; \
+  counters::_ctrs[i]._ctrNames = (char *)__FILE__;  \
   counters::_ctrs[i]._ctrLines = __LINE__; \
   salome_chrono aChrono##i(i);
 
@@ -67,7 +71,7 @@ protected:
 
 #else  // CHRONODEF
 
-#define CHRONO(i)
+#define CHRONO(i) using namespace SMDS;
 #define CHRONOSTOP(i)
 
 #endif // CHRONODEF
diff --git a/src/SMESHFiltersSelection/SMESH_NumberFilter.cxx b/src/SMESHFiltersSelection/SMESH_NumberFilter.cxx
index e5d480e35..78eed3a23 100644
--- a/src/SMESHFiltersSelection/SMESH_NumberFilter.cxx
+++ b/src/SMESHFiltersSelection/SMESH_NumberFilter.cxx
@@ -135,7 +135,7 @@ bool SMESH_NumberFilter::isOk (const SUIT_DataOwner* theDataOwner) const
   }
 
   // Verify number of sub-shapes
-  if (mySubShapeType == TopAbs_SHAPE);
+  if (mySubShapeType == TopAbs_SHAPE)
     return true;
 
   TopTools_IndexedMapOfShape aMap;
diff --git a/src/SMESHGUI/SMESHGUI.cxx b/src/SMESHGUI/SMESHGUI.cxx
index 27fe6571b..6ed14a435 100644
--- a/src/SMESHGUI/SMESHGUI.cxx
+++ b/src/SMESHGUI/SMESHGUI.cxx
@@ -4748,7 +4748,7 @@ void SMESHGUI::initialize( CAM_Application* app )
   popupMgr()->insert( separator(), -1, 0 );
   createPopupItem( SMESHOp::OpFileInformation,   OB, mesh, "&& selcount=1 && isImported" );
   createPopupItem( SMESHOp::OpMeshInformation,   OB, mesh_part );
-  createPopupItem( SMESHOp::OpFindElementByPoint,OB, mesh_group, "&& selcount=1" );
+  createPopupItem( SMESHOp::OpFindElementByPoint,OB, mesh_group, "&& selcount=1 && " + hasElems );
   createPopupItem( SMESHOp::OpOverallMeshQuality,OB, mesh_part );
   popupMgr()->insert( separator(), -1, 0 );
   createPopupItem( SMESHOp::OpCreateGroup,       OB, mesh, "&& selcount=1" );
@@ -4758,7 +4758,7 @@ void SMESHGUI::initialize( CAM_Application* app )
   createPopupItem( SMESHOp::OpEditHypothesis,    OB, hypo, "&& isEditableHyp");
   createPopupItem( SMESHOp::OpUnassign,          OB, hyp_alg );
   popupMgr()->insert( separator(), -1, 0 );
-  createPopupItem( SMESHOp::OpConvertMeshToQuadratic, OB, mesh_submesh );
+  createPopupItem( SMESHOp::OpConvertMeshToQuadratic, OB, mesh_submesh, "&& " + hasElems );
   createPopupItem( SMESHOp::OpCreateBoundaryElements, OB, mesh_group, "&& selcount=1 && dim>=2");
   popupMgr()->insert( separator(), -1, 0 );
 
@@ -4812,7 +4812,7 @@ void SMESHGUI::initialize( CAM_Application* app )
   popupMgr()->insert( separator(), -1, 0 );
   createPopupItem( SMESHOp::OpMeshInformation,    View, mesh_part );
   createPopupItem( SMESHOp::OpOverallMeshQuality, View, mesh_part );
-  createPopupItem( SMESHOp::OpFindElementByPoint, View, mesh );
+  createPopupItem( SMESHOp::OpFindElementByPoint, View, mesh, "&& " + hasElems);
   popupMgr()->insert( separator(), -1, 0 );
 
   createPopupItem( SMESHOp::OpUpdate,           OB + " " + View, mesh_part );
@@ -5994,6 +5994,9 @@ LightApp_Operation* SMESHGUI::createOperation( const int id ) const
   case SMESHOp::OpComputeSubMesh:
     op = new SMESHGUI_ComputeOp();
     break;
+  case SMESHOp::OpShowErrors:
+    op = new SMESHGUI_ShowErrorsOp();
+    break;
   case SMESHOp::OpPreCompute:
     op = new SMESHGUI_PrecomputeOp();
     break;
diff --git a/src/SMESHGUI/SMESHGUI_Displayer.cxx b/src/SMESHGUI/SMESHGUI_Displayer.cxx
index bfd048dc8..adeb2d736 100644
--- a/src/SMESHGUI/SMESHGUI_Displayer.cxx
+++ b/src/SMESHGUI/SMESHGUI_Displayer.cxx
@@ -89,7 +89,8 @@ SalomeApp_Study* SMESHGUI_Displayer::study() const
   return dynamic_cast<SalomeApp_Study*>( myApp->activeStudy() );
 }
 
-bool SMESHGUI_Displayer::canBeDisplayed( const QString& entry, const QString& viewer_type ) const {
+bool SMESHGUI_Displayer::canBeDisplayed( const QString& entry, const QString& viewer_type ) const
+{
   bool res = false;
   if(viewer_type != SVTK_Viewer::Type())
     return res;
@@ -97,25 +98,18 @@ bool SMESHGUI_Displayer::canBeDisplayed( const QString& entry, const QString& vi
   _PTR(SObject) obj = SMESH::getStudy()->FindObjectID( (const char*)entry.toUtf8() );
   CORBA::Object_var anObj = SMESH::SObjectToObject( obj );
   
-    /*
-    if( !CORBA::is_nil( anObj ) ) {
+  if ( !CORBA::is_nil( anObj ) )
+  {
     SMESH::SMESH_Mesh_var mesh = SMESH::SMESH_Mesh::_narrow( anObj );
-    if ( ! mesh->_is_nil() )
-    res = (mesh->NbNodes() > 0);
-    
-    SMESH::SMESH_subMesh_var aSubMeshObj = SMESH::SMESH_subMesh::_narrow( anObj );
-    if ( !aSubMeshObj->_is_nil() )
-    res = (aSubMeshObj->GetNumberOfNodes(true) > 0);
-    
-    SMESH::SMESH_GroupBase_var aGroupObj = SMESH::SMESH_GroupBase::_narrow( anObj );
-    if ( !aGroupObj->_is_nil() )
-    res = !aGroupObj->IsEmpty();
-    }*/
-  if( !CORBA::is_nil( anObj ) ) {
-    if(!SMESH::SMESH_Mesh::_narrow( anObj )->_is_nil() ||
-       !SMESH::SMESH_subMesh::_narrow( anObj )->_is_nil() ||
-       !SMESH::SMESH_GroupBase::_narrow( anObj )->_is_nil())
-    res = true;
+    if ( ! ( res = !mesh->_is_nil() ))
+    {
+      SMESH::SMESH_subMesh_var aSubMeshObj = SMESH::SMESH_subMesh::_narrow( anObj );
+      if ( ! ( res = !aSubMeshObj->_is_nil() ))
+      {
+        SMESH::SMESH_GroupBase_var aGroupObj = SMESH::SMESH_GroupBase::_narrow( anObj );
+        res = !aGroupObj->_is_nil();
+      }
+    }
   }
   return res;
 }
diff --git a/src/SMESHGUI/SMESHGUI_GroupDlg.cxx b/src/SMESHGUI/SMESHGUI_GroupDlg.cxx
index 5d33b4793..2313c768c 100644
--- a/src/SMESHGUI/SMESHGUI_GroupDlg.cxx
+++ b/src/SMESHGUI/SMESHGUI_GroupDlg.cxx
@@ -74,6 +74,7 @@
 #include <vtkActorCollection.h>
 
 // OCCT includes
+#include <TColStd_MapOfAsciiString.hxx>
 #include <TColStd_MapOfInteger.hxx>
 
 // Qt includes
@@ -522,19 +523,38 @@ QString SMESHGUI_GroupDlg::GetDefaultName(const QString& theOperation)
   return aName;
 }
 
+//=======================================================================
+//function : setDefaultName
+//purpose  : 
+//=======================================================================
+
 void  SMESHGUI_GroupDlg::setDefaultName() const
 {
   QString aResName;
-  int i=1;
-  QString aPrefix ="Group_";
-  _PTR(SObject) anObj;
-  do
+  const QString aPrefix ="Group_";
+
+  if ( myMesh->_is_nil() )
   {
-    aResName = aPrefix + QString::number( i++ );
-    anObj = SMESH::getStudy()->FindObject( aResName.toUtf8().constData() );
+    aResName = aPrefix + "1";
   }
-  while ( anObj );
-  myName->setText(aResName); 
+  else
+  {
+    SMESH::ListOfGroups_var allGroups = myMesh->GetGroups();
+    TColStd_MapOfAsciiString allNames( allGroups->length() );
+    for ( CORBA::ULong i = 0; i < allGroups->length(); ++i )
+    {
+      CORBA::String_var name = allGroups[i]->GetName();
+      allNames.Add( name.in() );
+    }
+    int i = 1;
+    while ( true )
+    {
+      aResName = aPrefix + QString::number( i++ );
+      if ( !allNames.Contains( aResName.toUtf8().constData() ))
+        break;
+    }
+  }
+  myName->setText(aResName);
 }
 
 //=================================================================================
@@ -944,7 +964,8 @@ bool SMESHGUI_GroupDlg::onApply()
 
     if (myGroup->_is_nil()) { // creation or conversion
       // check if group on geometry is not null
-      if (!myGroupOnGeom->_is_nil() || !myGroupOnFilter->_is_nil()) {
+      if (!myGroupOnGeom->_is_nil() || !myGroupOnFilter->_is_nil())
+      {
         if (myMesh->_is_nil())
           return false;
         if ( myGroupOnGeom->_is_nil() )
@@ -958,7 +979,8 @@ bool SMESHGUI_GroupDlg::onApply()
       }
     }
 
-    if (myGroup->_is_nil()) { // creation
+    if (myGroup->_is_nil()) // creation
+    {
       if (myMesh->_is_nil())
         return false;
 
@@ -990,10 +1012,10 @@ bool SMESHGUI_GroupDlg::onApply()
           myGroup->Add(anIdList.inout());
         }
       }
+    }
 
-
-    } else { // edition
-
+    else // edition
+    {
       resultGroup = SMESH::SMESH_GroupBase::_narrow( myGroup );
       isCreation = false;
 
@@ -1111,6 +1133,7 @@ bool SMESHGUI_GroupDlg::onApply()
     }
     anIsOk = true;
   }
+
   if (myGrpTypeId == 2) // group on filter
   {
     if ( myFilter->_is_nil() ) return false;
@@ -1170,7 +1193,7 @@ bool SMESHGUI_GroupDlg::onApply()
           Handle(SALOME_InteractiveObject) anIO = anActor->getIO();
           if ( isConversion ) { // need to reset TVisualObj and actor
             SMESH::RemoveVisualObjectWithActors( anIO->getEntry(), true );
-            SMESH::Update( anIO,true);
+            SMESH::Update( anIO, true );
             myActorsList.clear();
             anActor = SMESH::FindActorByEntry( anIO->getEntry() );
             if ( !anActor ) return false;
@@ -1402,9 +1425,11 @@ void SMESHGUI_GroupDlg::onObjectSelectionChanged()
       SALOME_ListIteratorOfListIO anIt (aList);
       for (; anIt.More(); anIt.Next())
       {
-        aGeomGroup = GEOMBase::ConvertIOinGEOMObject(anIt.Value());
-
+        CORBA::Object_var aGroupObj = SMESH::IObjectToObject(anIt.Value());
+        if (CORBA::is_nil(aGroupObj))
+          continue;
         // Check if the object is a geometry group
+        aGeomGroup = GEOM::GEOM_Object::_narrow(aGroupObj);
         if (CORBA::is_nil(aGeomGroup))
           continue;
 
diff --git a/src/SMESHGUI/SMESHGUI_Selection.cxx b/src/SMESHGUI/SMESHGUI_Selection.cxx
index fdf9a8d1c..77d5abaa2 100644
--- a/src/SMESHGUI/SMESHGUI_Selection.cxx
+++ b/src/SMESHGUI/SMESHGUI_Selection.cxx
@@ -43,8 +43,6 @@
 #include <SalomeApp_Application.h>
 #include <LightApp_VTKSelector.h>
 #include <SVTK_ViewWindow.h>
-#include <LightApp_SelectionMgr.h>
-#include <LightApp_DataOwner.h>
 
 // IDL includes
 #include <SALOMEconfig.h>
@@ -79,8 +77,6 @@ void SMESHGUI_Selection::init( const QString& client, LightApp_SelectionMgr* mgr
 
   if( mgr )
   {
-    myOwners.clear();
-    mgr->selected(myOwners, client);
     for( int i=0, n=count(); i<n; i++ ) {
       myTypes.append( typeName( type( entry( i ) ) ) );
       myControls.append( controlMode( i ) );
@@ -659,20 +655,22 @@ bool SMESHGUI_Selection::hasGeomReference( int ind ) const
 
 bool SMESHGUI_Selection::canBreakLink( int ind ) const
 {
-  if ( ind >= 0 && ind < myTypes.count()) {
-    if (isReference(ind)) {
-      SUIT_DataOwner* aOwn = myOwners.at(ind);
-      LightApp_DataOwner* sowner = dynamic_cast<LightApp_DataOwner*>(aOwn);
-      QString aEntry = sowner->entry();
-      _PTR(SObject) aSObject = SMESH::getStudy()->FindObjectID(aEntry.toStdString());
-      _PTR(SObject) aFatherObj = aSObject->GetFather();
-      _PTR(SComponent) aComponent = aFatherObj->GetFatherComponent();
-      if (aComponent->ComponentDataType() == "SMESH") {
-        QString aObjEntry = entry(ind);
-        _PTR(SObject) aGeomSObject = SMESH::getStudy()->FindObjectID(aObjEntry.toStdString());
-        GEOM::GEOM_Object_var aObject = SMESH::SObjectToInterface<GEOM::GEOM_Object>(aGeomSObject);
-        if (!aObject->_is_nil())
-          return aObject->IsParametrical();
+  if ( ind >= 0 && isReference(ind) ) {
+    QString aEntry = objectInfo(ind, OI_RefEntry).toString();
+    if (!aEntry.isEmpty()) {
+      _PTR(SObject) aSObject = SMESH::getStudy()->FindObjectID( aEntry.toStdString());
+      if (aSObject) {
+        _PTR(SObject) aFatherObj = aSObject->GetFather();
+        if (aFatherObj) {
+          _PTR(SComponent) aComponent = aFatherObj->GetFatherComponent();
+          if (aComponent && aComponent->ComponentDataType() == "SMESH") {
+            QString aObjEntry = entry(ind);
+            _PTR(SObject) aGeomSObject = SMESH::getStudy()->FindObjectID(aObjEntry.toStdString());
+            GEOM::GEOM_Object_var aObject = SMESH::SObjectToInterface<GEOM::GEOM_Object>(aGeomSObject);
+            if (!aObject->_is_nil())
+              return aObject->IsParametrical();
+          }
+        }
       }
     }
   }
diff --git a/src/SMESHGUI/SMESHGUI_Selection.h b/src/SMESHGUI/SMESHGUI_Selection.h
index b9727108f..710cfcf83 100644
--- a/src/SMESHGUI/SMESHGUI_Selection.h
+++ b/src/SMESHGUI/SMESHGUI_Selection.h
@@ -32,7 +32,6 @@
 
 // SALOME GUI includes
 #include <LightApp_Selection.h>
-#include <SUIT_DataOwner.h>
 
 // SALOME KERNEL includes
 #include <SALOMEDSClient_definitions.hxx>
@@ -99,7 +98,6 @@ private:
   QStringList             myTypes;
   QStringList             myControls;
   QList<SMESH_Actor*>     myActors;
-  SUIT_DataOwnerPtrList   myOwners;
 };
 
 #endif // SMESHGUI_SELECTION_H
diff --git a/src/SMESH_I/SMESH_2smeshpy.cxx b/src/SMESH_I/SMESH_2smeshpy.cxx
index 9ad385290..b36b31c13 100644
--- a/src/SMESH_I/SMESH_2smeshpy.cxx
+++ b/src/SMESH_I/SMESH_2smeshpy.cxx
@@ -3939,7 +3939,7 @@ int _pyCommand::GetArgBeginning() const
   if ( pos == UNKNOWN )
   {
     pos = GetBegPos( METHOD_IND ) + myMeth.Length();
-    if ( pos < 1 )
+    if ( pos < 1 && Length() >= 4 )
       pos = myString.Location( "(", 4, Length() ); // 4 = strlen("b.c(")
   }
   return pos;
@@ -4024,6 +4024,11 @@ TCollection_AsciiString _pyCommand::GetWord( const _AString & theString,
   }
   theStartPos = beg;
   //cout << theString << " ---- " << beg << " - " << end << endl;
+  if ( end > theString.Length() )
+  {
+    theStartPos = EMPTY;
+    return theEmptyString;
+  }
   return theString.SubString( beg, end );
 }
 
diff --git a/src/SMESH_I/SMESH_Gen_i.hxx b/src/SMESH_I/SMESH_Gen_i.hxx
index 895f745fc..d6cb5b130 100644
--- a/src/SMESH_I/SMESH_Gen_i.hxx
+++ b/src/SMESH_I/SMESH_Gen_i.hxx
@@ -494,6 +494,8 @@ public:
 
   void CleanPythonTrace();
 
+  static int CountInPyDump(const TCollection_AsciiString& text);
+
   // *****************************************
   // Internal methods
   // *****************************************
diff --git a/src/SMESH_I/SMESH_Mesh_i.cxx b/src/SMESH_I/SMESH_Mesh_i.cxx
index cf2a16ea6..4c2ca27c3 100644
--- a/src/SMESH_I/SMESH_Mesh_i.cxx
+++ b/src/SMESH_I/SMESH_Mesh_i.cxx
@@ -2058,12 +2058,12 @@ void SMESH_Mesh_i::ReplaceShape(GEOM::GEOM_Object_ptr theNewGeom)
   TopoDS_Shape S = _impl->GetShapeToMesh();
   GEOM_Client* geomClient = _gen_i->GetShapeReader();
   TCollection_AsciiString aIOR;
-  CORBA::String_var ior;
   if ( geomClient->Find( S, aIOR ))
     geomClient->RemoveShapeFromBuffer( aIOR );
 
   // clear buffer also for sub-groups
-  const std::set<SMESHDS_GroupBase*>& groups = _impl->GetMeshDS()->GetGroups();
+  SMESHDS_Mesh* meshDS = _impl->GetMeshDS();
+  const std::set<SMESHDS_GroupBase*>& groups = meshDS->GetGroups();
   std::set<SMESHDS_GroupBase*>::const_iterator g = groups.begin();
   for (; g != groups.end(); ++g)
     if (const SMESHDS_GroupOnGeom* group = dynamic_cast<SMESHDS_GroupOnGeom*>(*g))
@@ -2073,12 +2073,21 @@ void SMESH_Mesh_i::ReplaceShape(GEOM::GEOM_Object_ptr theNewGeom)
         geomClient->RemoveShapeFromBuffer( aIOR );
     }
 
+  // clear buffer also for sub-meshes
+  std::map<int, SMESH_subMesh_i*>::const_iterator aSubMeshIter = _mapSubMesh_i.cbegin();
+  for(; aSubMeshIter != _mapSubMesh_i.cend(); aSubMeshIter++) {
+    int aShapeID = aSubMeshIter->first;
+    const TopoDS_Shape& aSubShape = meshDS->IndexToShape(aShapeID);
+    TCollection_AsciiString aShapeIOR;
+    if ( geomClient->Find( aSubShape, aShapeIOR ))
+      geomClient->RemoveShapeFromBuffer( aShapeIOR );
+  }
+
   typedef struct {
     int shapeID, fromID, toID; // indices of elements of a sub-mesh
   } TRange;
   std::vector< TRange > elemRanges, nodeRanges; // elements of sub-meshes
   std::vector< SMDS_PositionPtr > positions; // node positions
-  SMESHDS_Mesh* meshDS = _impl->GetMeshDS();
   if ( !geomChanged )
   {
     // store positions of elements on geometry
@@ -5821,15 +5830,21 @@ void SMESH_Mesh_i::CreateGroupServants()
     GEOM::GEOM_Object_var shapeVar = _gen_i->ShapeToGeomObject( shape );
     _gen_i->PublishGroup( aMesh, groupVar, shapeVar, group->GetName());
   }
+
   if ( !addedIDs.empty() )
   {
     // python dump
-    set<int>::iterator id = addedIDs.begin();
-    for ( ; id != addedIDs.end(); ++id )
+    map<int, SMESH::SMESH_GroupBase_ptr>::iterator i_grp = _mapGroups.begin();
+    for ( int index = 0; i_grp != _mapGroups.end(); ++index, ++i_grp )
     {
-      map<int, SMESH::SMESH_GroupBase_ptr>::iterator it = _mapGroups.find(*id);
-      int i = std::distance( _mapGroups.begin(), it );
-      TPythonDump() << it->second << " = " << aMesh << ".GetGroups()[ "<< i << " ]";
+      set<int>::iterator it = addedIDs.find( i_grp->first );
+      if ( it != addedIDs.end() )
+      {
+        TPythonDump() << i_grp->second << " = " << aMesh << ".GetGroups()[ "<< index << " ]";
+        addedIDs.erase( it );
+        if ( addedIDs.empty() )
+          break;
+      }
     }
   }
 }
@@ -6757,7 +6772,8 @@ TListOfListOfInt SMESH_Mesh_i::findConcurrentSubMeshes()
   TDimHypList dimHypListArr[4]; // dimHyp list for each shape dimension
 
   map<int, ::SMESH_subMesh*>::iterator i_sm = _mapSubMesh.begin();
-  for ( ; i_sm != _mapSubMesh.end(); i_sm++ ) {
+  for ( ; i_sm != _mapSubMesh.end(); i_sm++ )
+  {
     ::SMESH_subMesh* sm = (*i_sm).second;
     // shape of submesh
     const TopoDS_Shape& aSubMeshShape = sm->GetSubShape();
@@ -6767,7 +6783,8 @@ TListOfListOfInt SMESH_Mesh_i::findConcurrentSubMeshes()
     // Find out dimensions where the submesh can be concurrent.
     // We define the dimensions by algo of each of hypotheses in hypList
     list <const SMESHDS_Hypothesis*>::const_iterator hypIt = hypList.begin();
-    for( ; hypIt != hypList.end(); hypIt++ ) {
+    for( ; hypIt != hypList.end(); hypIt++ )
+    {
       SMESH_Algo* anAlgo = 0;
       const SMESH_Hypothesis* hyp = dynamic_cast<const SMESH_Hypothesis*>(*hypIt);
       if ( hyp->GetType() != SMESHDS_Hypothesis::PARAM_ALGO )
diff --git a/src/SMESH_I/SMESH_PythonDump.cxx b/src/SMESH_I/SMESH_PythonDump.cxx
index 04a0fc9a8..a09766949 100644
--- a/src/SMESH_I/SMESH_PythonDump.cxx
+++ b/src/SMESH_I/SMESH_PythonDump.cxx
@@ -40,6 +40,8 @@
 #include <TColStd_HSequenceOfInteger.hxx>
 #include <TCollection_AsciiString.hxx>
 
+#include <cstring>
+
 #ifdef _DEBUG_
 static int MYDEBUG = 0;
 #else
@@ -77,7 +79,7 @@ namespace SMESH
       TCollection_AsciiString aCollection(Standard_CString(aString.c_str()));
       if(!aCollection.IsEmpty())
       {
-        const std::string & objEntry = SMESH_Gen_i::GetSMESHGen()->GetLastObjEntry();
+        const std::string & objEntry = aSMESHGen->GetLastObjEntry();
         if ( !objEntry.empty() )
           aCollection += (TVar::ObjPrefix() + objEntry ).c_str();
         aSMESHGen->AddToPythonScript(aCollection);
@@ -1078,10 +1080,15 @@ TCollection_AsciiString SMESH_Gen_i::DumpPython_impl
   TCollection_AsciiString anUpdatedScript;
 
   Resource_DataMapOfAsciiStringAsciiString mapRemoved;
-  Resource_DataMapOfAsciiStringAsciiString mapEntries; // names and entries present in anUpdatedScript
+  Resource_DataMapOfAsciiStringAsciiString mapEntries; // { entry: name } present in anUpdatedScript
   Standard_Integer objectCounter = 0;
   TCollection_AsciiString anEntry, aName, aGUIName, aBaseName("smeshObj_");
 
+  std::string           compDataType = ComponentDataType(); // SMESH module's data type
+  SALOMEDS::SComponent_var   smeshSO = getStudyServant()->FindComponent( compDataType.c_str() );
+  CORBA::String_var          smeshID = smeshSO->GetID();
+  TCollection_AsciiString smeshEntry = smeshID.in();
+
   // Treat every script line and add it to anUpdatedScript
   for ( linesIt = lines.begin(); linesIt != lines.end(); ++linesIt )
   {
@@ -1133,15 +1140,22 @@ TCollection_AsciiString SMESH_Gen_i::DumpPython_impl
         }
         else
         {
-          // Removed Object
-          do {
-            aName = aBaseName + (++objectCounter);
-          } while (theObjectNames.IsBound(aName));
+          if ( !anEntry.StartsWith( smeshEntry )) // not SMESH object
+          {
+            aName = SMESH::TPythonDump::NotPublishedObjectName();
+          }
+          else
+          {
+            // Removed Object
+            do {
+              aName = aBaseName + (++objectCounter);
+            } while (theObjectNames.IsBound(aName));
 
-          if ( !aRemovedObjIDs.count( anEntry ) && aLine.Value(1) != '#')
-            mapRemoved.Bind(anEntry, aName);
+            if ( !aRemovedObjIDs.count( anEntry ) && aLine.Value(1) != '#')
+              mapRemoved.Bind(anEntry, aName);
 
-          theObjectNames.Bind(anEntry, aName);
+            theObjectNames.Bind(anEntry, aName);
+          }
         }
         theObjectNames.Bind(aName, anEntry); // to detect same name of diff objects
       }
@@ -1388,3 +1402,55 @@ void SMESH_Gen_i::CleanPythonTrace()
     myPythonScript->Clear();
   }
 }
+
+//================================================================================
+/*!
+ * \brief Count inclusions of a string in a raw Python dump script
+ */
+//================================================================================
+
+int SMESH_Gen_i::CountInPyDump(const TCollection_AsciiString& theText)
+{
+  int count = 0;
+
+  SALOMEDS::Study_var aStudy = getStudyServant();
+  if ( CORBA::is_nil( aStudy ))
+    return count;
+
+  SMESH_Gen_i* me = GetSMESHGen();
+  CORBA::String_var compDataType = me->ComponentDataType();
+  SALOMEDS::SObject_wrap aSO = aStudy->FindComponent( compDataType.in() );
+  if ( CORBA::is_nil( aSO ))
+    return count;
+
+  // Trace saved in the study
+  SALOMEDS::GenericAttribute_wrap attr;
+  if ( aSO->FindAttribute( attr.inout(), "AttributePythonObject" ))
+  {
+    SALOMEDS::AttributePythonObject_var pyAttr =
+      SALOMEDS::AttributePythonObject::_narrow( attr );
+    CORBA::String_var script = pyAttr->GetObject();
+    for ( const char * scriptPos = script.in(); true; ++scriptPos )
+      if (( scriptPos = strstr( scriptPos, theText.ToCString() )))
+        ++count;
+      else
+        break;
+  }
+
+  // New python commands
+  if ( !me->myPythonScript.IsNull() )
+  {
+    const int nbLines = me->myPythonScript->Length();
+    for ( int i = 1; i <= nbLines; ++i )
+    {
+      const TCollection_AsciiString& line = me->myPythonScript->Value( i );
+      for ( int loc = 1; loc <= line.Length(); ++loc )
+        if (( loc = line.Location( theText, loc, line.Length() )))
+          ++count;
+        else
+          break;
+    }
+  }
+
+  return count;
+}
diff --git a/src/SMESH_SWIG/smeshBuilder.py b/src/SMESH_SWIG/smeshBuilder.py
index cb0191fcc..3cb09bd0e 100644
--- a/src/SMESH_SWIG/smeshBuilder.py
+++ b/src/SMESH_SWIG/smeshBuilder.py
@@ -2636,7 +2636,15 @@ class Mesh(metaclass = MeshMeta):
         elif tgeo == "SOLID" or tgeo == "COMPSOLID":
             typ = VOLUME
         elif tgeo == "COMPOUND":
-            sub = self.geompyD.SubShapeAll( shape, self.geompyD.ShapeType["SHAPE"])
+            try:
+              sub = self.geompyD.SubShapeAll( shape, self.geompyD.ShapeType["SHAPE"])
+            except:
+              # try to get the SHAPERSTUDY engine directly, because GetGen does not work because of
+              # simplification of access in geomBuilder: omniORB.registerObjref
+              from SHAPERSTUDY_utils import getEngine
+              gen = getEngine()
+              if gen:
+                sub = gen.GetIShapesOperations().ExtractSubShapes(shape, self.geompyD.ShapeType["SHAPE"], False)
             if not sub:
                 raise ValueError("_groupTypeFromShape(): empty geometric group or compound '%s'" % GetName(shape))
             return self._groupTypeFromShape( sub[0] )
diff --git a/src/SMESH_SWIG/smesh_algorithm.py b/src/SMESH_SWIG/smesh_algorithm.py
index e7003b509..8ee3718a3 100644
--- a/src/SMESH_SWIG/smesh_algorithm.py
+++ b/src/SMESH_SWIG/smesh_algorithm.py
@@ -340,9 +340,20 @@ class Mesh_Algorithm:
         if faces and isinstance( faces[0], geomBuilder.GEOM._objref_GEOM_Object ):
             faceIDs = []
             for shape in faces:
-                ff = self.mesh.geompyD.SubShapeAll( shape, self.mesh.geompyD.ShapeType["FACE"] )
-                for f in ff:
+                try:
+                  ff = self.mesh.geompyD.SubShapeAll( shape, self.mesh.geompyD.ShapeType["FACE"] )
+                  for f in ff:
                     faceIDs.append( self.mesh.geompyD.GetSubShapeID(self.mesh.geom, f))
+                except:
+                  # try to get the SHAPERSTUDY engine directly, because GetGen does not work because of
+                  # simplification of access in geomBuilder: omniORB.registerObjref
+                  from SHAPERSTUDY_utils import getEngine
+                  gen = getEngine()
+                  if gen:
+                    aShapeOp = gen.GetIShapesOperations()
+                    ff = aShapeOp.ExtractSubShapes( shape, self.mesh.geompyD.ShapeType["FACE"], False)
+                    for f in ff:
+                      faceIDs.append( aShapeOp.GetSubShapeIndex( self.mesh.geom, f ))
             faces = faceIDs
         hyp = self.Hypothesis("ViscousLayers",
                               [thickness, numberOfLayers, stretchFactor, faces, isFacesToIgnore],
@@ -392,9 +403,20 @@ class Mesh_Algorithm:
         if edges and isinstance( edges[0], geomBuilder.GEOM._objref_GEOM_Object ):
             edgeIDs = []
             for shape in edges:
+              try:
                 ee = self.mesh.geompyD.SubShapeAll( shape, self.mesh.geompyD.ShapeType["EDGE"])
                 for e in ee:
-                    edgeIDs.append( self.mesh.geompyD.GetSubShapeID( self.mesh.geom, e ))
+                  edgeIDs.append( self.mesh.geompyD.GetSubShapeID( self.mesh.geom, e ))
+              except:
+                # try to get the SHAPERSTUDY engine directly, because GetGen does not work because of
+                # simplification of access in geomBuilder: omniORB.registerObjref
+                from SHAPERSTUDY_utils import getEngine
+                gen = getEngine()
+                if gen:
+                  aShapeOp = gen.GetIShapesOperations()
+                  ee = aShapeOp.ExtractSubShapes( shape, self.mesh.geompyD.ShapeType["EDGE"], False)
+                  for e in ee:
+                    edgeIDs.append( aShapeOp.GetSubShapeIndex( self.mesh.geom, e ))
             edges = edgeIDs
         hyp = self.Hypothesis("ViscousLayers2D",
                               [thickness, numberOfLayers, stretchFactor, edges, isEdgesToIgnore],
diff --git a/src/StdMeshers/StdMeshers_Cartesian_3D.cxx b/src/StdMeshers/StdMeshers_Cartesian_3D.cxx
index ab129db36..6f284e4ba 100644
--- a/src/StdMeshers/StdMeshers_Cartesian_3D.cxx
+++ b/src/StdMeshers/StdMeshers_Cartesian_3D.cxx
@@ -362,6 +362,9 @@ namespace
     gp_XYZ             _origin;
     gp_Mat             _invB; // inverted basis of _axes
 
+    // index shift within _nodes of nodes of a cell from the 1st node
+    int                _nodeShift[8];
+
     vector< const SMDS_MeshNode* >    _nodes; // mesh nodes at grid nodes
     vector< const F_IntersectPoint* > _gridIntP; // grid node intersection with geometry
     ObjectPool< E_IntersectPoint >    _edgeIntPool; // intersections with EDGEs
@@ -427,6 +430,7 @@ namespace
     bool IsBoundaryFace( TGeomID face ) const { return _geometry._boundaryFaces.Contains( face ); }
     void SetOnShape( const SMDS_MeshNode* n, const F_IntersectPoint& ip, bool unset=false );
     bool IsToCheckNodePos() const { return !_toAddEdges && _toCreateFaces; }
+    bool IsToRemoveExcessEntities() const { return !_toAddEdges; }
 
     void SetCoordinates(const vector<double>& xCoords,
                         const vector<double>& yCoords,
@@ -442,12 +446,14 @@ namespace
    */
   struct CellsAroundLink
   {
+    int    _iDir;
     int    _dInd[4][3];
     size_t _nbCells[3];
     int    _i,_j,_k;
     Grid*  _grid;
 
     CellsAroundLink( Grid* grid, int iDir ):
+      _iDir( iDir ),
       _dInd{ {0,0,0}, {0,0,0}, {0,0,0}, {0,0,0} },
       _nbCells{ grid->_coords[0].size() - 1,
           grid->_coords[1].size() - 1,
@@ -466,7 +472,7 @@ namespace
       _j = j - _dInd[iL][1];
       _k = k - _dInd[iL][2];
     }
-    bool GetCell( int iL, int& i, int& j, int& k, int& cellIndex )
+    bool GetCell( int iL, int& i, int& j, int& k, int& cellIndex, int& linkIndex )
     {
       i =  _i + _dInd[iL][0];
       j =  _j + _dInd[iL][1];
@@ -476,6 +482,7 @@ namespace
            k < 0 || k >= (int)_nbCells[2] )
         return false;
       cellIndex = _grid->CellIndex( i,j,k );
+      linkIndex = iL + _iDir * 4;
       return true;
     }
   };
@@ -724,11 +731,11 @@ namespace
     {
       struct _Split // data of a link split
       {
-        int    _linkID;      // hex link ID
+        int    _linkID;          // hex link ID
         _Node* _nodes[2];
         int    _iCheckIteration; // iteration where split is tried as Hexahedron split
         _Link* _checkedSplit;    // split set to hex links
-        bool   _isUsed;    // used in a volume
+        bool   _isUsed;          // used in a volume
 
         _Split( _Link & split, int iLink ):
           _linkID( iLink ), _nodes{ split._nodes[0], split._nodes[1] },
@@ -759,9 +766,13 @@ namespace
     // --------------------------------------------------------------------------------
     struct _Face
     {
+      SMESH_Block::TShapeID   _name;
       vector< _OrientedLink > _links;       // links on GridLine's
       vector< _Link >         _polyLinks;   // links added to close a polygonal face
       vector< _Node* >        _eIntNodes;   // nodes at intersection with EDGEs
+
+      _Face():_name( SMESH_Block::ID_NONE )
+      {}
       bool IsPolyLink( const _OrientedLink& ol )
       {
         return _polyLinks.empty() ? false :
@@ -789,41 +800,79 @@ namespace
     // --------------------------------------------------------------------------------
     struct _volumeDef // holder of nodes of a volume mesh element
     {
+      typedef void* _ptr;
+
       struct _nodeDef
       {
         const SMDS_MeshNode*    _node; // mesh node at hexahedron corner
         const B_IntersectPoint* _intPoint;
 
+        _nodeDef(): _node(0), _intPoint(0) {}
         _nodeDef( _Node* n ): _node( n->_node), _intPoint( n->_intPoint ) {}
         const SMDS_MeshNode*    Node() const
         { return ( _intPoint && _intPoint->_node ) ? _intPoint->_node : _node; }
         const E_IntersectPoint* EdgeIntPnt() const
         { return static_cast< const E_IntersectPoint* >( _intPoint ); }
+        _ptr Ptr() const { return Node() ? (_ptr) Node() : (_ptr) EdgeIntPnt(); }
+        bool operator==(const _nodeDef& other ) const { return Ptr() == other.Ptr(); }
       };
+
       vector< _nodeDef >      _nodes;
       vector< int >           _quantities;
       _volumeDef*             _next; // to store several _volumeDefs in a chain
       TGeomID                 _solidID;
       const SMDS_MeshElement* _volume; // new volume
 
+      vector< SMESH_Block::TShapeID > _names; // name of side a polygon originates from
+
       _volumeDef(): _next(0), _solidID(0), _volume(0) {}
       ~_volumeDef() { delete _next; }
       _volumeDef( _volumeDef& other ):
         _next(0), _solidID( other._solidID ), _volume( other._volume )
-      { _nodes.swap( other._nodes ); _quantities.swap( other._quantities ); other._volume = 0; }
+      { _nodes.swap( other._nodes ); _quantities.swap( other._quantities ); other._volume = 0;
+        _names.swap( other._names ); }
 
-      void Set( const vector< _Node* >& nodes, const vector< int >& quant = vector< int >() )
-      { _nodes.assign( nodes.begin(), nodes.end() ); _quantities = quant; }
+      size_t size() const { return 1 + ( _next ? _next->size() : 0 ); }
+      _volumeDef* at(int index)
+      { return index == 0 ? this : ( _next ? _next->at(index-1) : _next ); }
 
       void Set( _Node** nodes, int nb )
       { _nodes.assign( nodes, nodes + nb ); }
 
       void SetNext( _volumeDef* vd )
       { if ( _next ) { _next->SetNext( vd ); } else { _next = vd; }}
+
+      bool IsEmpty() const { return (( _nodes.empty() ) &&
+                                     ( !_next || _next->IsEmpty() )); }
+      bool IsPolyhedron() const { return ( !_quantities.empty() ||
+                                           ( _next && !_next->_quantities.empty() )); }
+
+
+      struct _linkDef: public std::pair<_ptr,_ptr> // to join polygons in removeExcessSideDivision()
+      {
+        _nodeDef _node1;//, _node2;
+        mutable /*const */_linkDef *_prev, *_next;
+        size_t _loopIndex;
+
+        _linkDef():_prev(0), _next(0) {}
+
+        void init( const _nodeDef& n1, const _nodeDef& n2, size_t iLoop )
+        {
+          _node1     = n1; //_node2 = n2;
+          _loopIndex = iLoop;
+          first      = n1.Ptr();
+          second     = n2.Ptr();
+          if ( first > second ) std::swap( first, second );
+        }
+        void setNext( _linkDef* next )
+        {
+          _next = next;
+          next->_prev = this;
+        }
+      };
     };
 
     // topology of a hexahedron
-    int   _nodeShift[8];
     _Node _hexNodes [8];
     _Link _hexLinks [12];
     _Face _hexQuads [6];
@@ -837,7 +886,7 @@ namespace
     // additional nodes created at intersection points
     vector< _Node > _intNodes;
 
-    // nodes inside the hexahedron (at VERTEXes)
+    // nodes inside the hexahedron (at VERTEXes) refer to _intNodes
     vector< _Node* > _vIntNodes;
 
     // computed volume elements
@@ -858,7 +907,7 @@ namespace
     Hexahedron(Grid* grid);
     int MakeElements(SMESH_MesherHelper&                      helper,
                      const map< TGeomID, vector< TGeomID > >& edge2faceIDsMap);
-    void ComputeElements( const Solid* solid = 0, int solidIndex = -1 );
+    void computeElements( const Solid* solid = 0, int solidIndex = -1 );
 
   private:
     Hexahedron(const Hexahedron& other, size_t i, size_t j, size_t k, int cellID );
@@ -893,6 +942,9 @@ namespace
                       const map< TGeomID, vector< TGeomID > >& edge2faceIDsMap );
     void getVolumes( vector< const SMDS_MeshElement* > & volumes );
     void getBoundaryElems( vector< const SMDS_MeshElement* > & boundaryVolumes );
+    void removeExcessSideDivision(const vector< Hexahedron* >& allHexa);
+    void removeExcessNodes(vector< Hexahedron* >& allHexa);
+    void preventVolumesOverlapping();
     TGeomID getAnyFace() const;
     void cutByExtendedInternal( std::vector< Hexahedron* >& hexes,
                                 const TColStd_MapOfInteger& intEdgeIDs );
@@ -927,7 +979,7 @@ namespace
       TID2Nb::iterator id2nb = id2nbMap.insert( s0 ).first;
       id2nb->second++;
     }
-  };
+  }; // class Hexahedron
 
 #ifdef WITH_TBB
   // --------------------------------------------------------------------------
@@ -942,7 +994,7 @@ namespace
     {
       for ( size_t i = r.begin(); i != r.end(); ++i )
         if ( Hexahedron* hex = _hexVec[ i ] )
-          hex->ComputeElements();
+          hex->computeElements();
     }
   };
   // --------------------------------------------------------------------------
@@ -2089,14 +2141,14 @@ namespace
     size_t i101 = i100 + dz;
     size_t i011 = i010 + dz;
     size_t i111 = i110 + dz;
-    _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V000 )] = i000;
-    _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V100 )] = i100;
-    _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V010 )] = i010;
-    _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V110 )] = i110;
-    _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V001 )] = i001;
-    _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V101 )] = i101;
-    _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V011 )] = i011;
-    _nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V111 )] = i111;
+    grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V000 )] = i000;
+    grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V100 )] = i100;
+    grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V010 )] = i010;
+    grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V110 )] = i110;
+    grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V001 )] = i001;
+    grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V101 )] = i101;
+    grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V011 )] = i011;
+    grid->_nodeShift[ SMESH_Block::ShapeIndex( SMESH_Block::ID_V111 )] = i111;
 
     vector< int > idVec;
     // set nodes to links
@@ -2112,8 +2164,10 @@ namespace
     int interlace[4] = { 0, 3, 1, 2 }; // to walk by links around a face: { u0, 1v, u1, 0v }
     for ( int faceID = SMESH_Block::ID_Fxy0; faceID <= SMESH_Block::ID_F1yz; ++faceID )
     {
-      SMESH_Block::GetFaceEdgesIDs( faceID, idVec );
       _Face& quad = _hexQuads[ SMESH_Block::ShapeIndex( faceID )];
+      quad._name = (SMESH_Block::TShapeID) faceID;
+
+      SMESH_Block::GetFaceEdgesIDs( faceID, idVec );
       bool revFace = ( faceID == SMESH_Block::ID_Fxy0 ||
                        faceID == SMESH_Block::ID_Fx1z ||
                        faceID == SMESH_Block::ID_F0yz );
@@ -2141,9 +2195,6 @@ namespace
     _polygons.reserve(100); // to avoid reallocation;
 
     // copy topology
-    for ( int i = 0; i < 8; ++i )
-      _nodeShift[i] = other._nodeShift[i];
-
     for ( int i = 0; i < 12; ++i )
     {
       const _Link& srcLink = other._hexLinks[ i ];
@@ -2156,6 +2207,7 @@ namespace
     {
       const _Face& srcQuad = other._hexQuads[ i ];
       _Face&       tgtQuad = this->_hexQuads[ i ];
+      tgtQuad._name = srcQuad._name;
       tgtQuad._links.resize(4);
       for ( int j = 0; j < 4; ++j )
       {
@@ -2188,8 +2240,8 @@ namespace
     _origNodeInd = _grid->NodeIndex( _i,_j,_k );
     for ( int iN = 0; iN < 8; ++iN )
     {
-      _hexNodes[iN]._node     = _grid->_nodes   [ _origNodeInd + _nodeShift[iN] ];
-      _hexNodes[iN]._intPoint = _grid->_gridIntP[ _origNodeInd + _nodeShift[iN] ];
+      _hexNodes[iN]._node     = _grid->_nodes   [ _origNodeInd + _grid->_nodeShift[iN] ];
+      _hexNodes[iN]._intPoint = _grid->_gridIntP[ _origNodeInd + _grid->_nodeShift[iN] ];
 
       if ( _hexNodes[iN]._intPoint ) // intersection with a FACE
       {
@@ -2337,6 +2389,7 @@ namespace
   {
     _i = i; _j = j; _k = k;
 
+    bool isCompute = solid;
     if ( !solid )
       solid = _grid->GetSolid();
 
@@ -2348,8 +2401,8 @@ namespace
     {
       _hexNodes[iN]._isInternalFlags = 0;
 
-      _hexNodes[iN]._node     = _grid->_nodes   [ _origNodeInd + _nodeShift[iN] ];
-      _hexNodes[iN]._intPoint = _grid->_gridIntP[ _origNodeInd + _nodeShift[iN] ];
+      _hexNodes[iN]._node     = _grid->_nodes   [ _origNodeInd + _grid->_nodeShift[iN] ];
+      _hexNodes[iN]._intPoint = _grid->_gridIntP[ _origNodeInd + _grid->_nodeShift[iN] ];
 
       if ( _hexNodes[iN]._node && !solid->Contains( _hexNodes[iN]._node->GetShapeID() ))
         _hexNodes[iN]._node = 0;
@@ -2366,6 +2419,9 @@ namespace
     _intNodes.clear();
     _vIntNodes.clear();
 
+    if ( !isCompute )
+      return;
+
     if ( _nbFaceIntNodes + _eIntPoints.size()                  > 0 &&
          _nbFaceIntNodes + _eIntPoints.size() + _nbCornerNodes > 3)
     {
@@ -2545,7 +2601,8 @@ namespace
       } // loop on _eIntPoints
     }
 
-    else if ( 3 < _nbCornerNodes && _nbCornerNodes < 8 ) // _nbFaceIntNodes == 0
+    else if (( 3 < _nbCornerNodes && _nbCornerNodes < 8 ) || // _nbFaceIntNodes == 0
+             ( !_grid->_geometry.IsOneSolid() ))
     {
       _Link split;
       // create sub-links (_splits) of whole links
@@ -2569,7 +2626,7 @@ namespace
   /*!
    * \brief Compute mesh volumes resulted from intersection of the Hexahedron
    */
-  void Hexahedron::ComputeElements( const Solid* solid, int solidIndex )
+  void Hexahedron::computeElements( const Solid* solid, int solidIndex )
   {
     if ( !solid )
     {
@@ -2583,7 +2640,7 @@ namespace
           for ( size_t i = 0; i < nbSolids; ++i )
           {
             solid = _grid->GetSolid( solidIDs[i] );
-            ComputeElements( solid, i );
+            computeElements( solid, i );
             if ( !_volumeDefs._nodes.empty() && i < nbSolids - 1 )
               _volumeDefs.SetNext( new _volumeDef( _volumeDefs ));
           }
@@ -2632,6 +2689,9 @@ namespace
     for ( int iN = 0; iN < 8; ++iN )
       _hexNodes[iN]._usedInFace = 0;
 
+    if ( intFlag & IS_CUT_BY_INTERNAL_FACE && !_grid->_toAddEdges ) // Issue #19913
+      preventVolumesOverlapping();
+
     // Create polygons from quadrangles
     // --------------------------------
 
@@ -2639,7 +2699,8 @@ namespace
     vector<_Node*>          chainNodes;
     _Face*                  coplanarPolyg;
 
-    bool hasEdgeIntersections = !_eIntPoints.empty();
+    const bool hasEdgeIntersections = !_eIntPoints.empty();
+    const bool toCheckSideDivision = isImplementEdges() || intFlag & IS_CUT_BY_INTERNAL_FACE;
 
     for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides of a hexahedron
     {
@@ -2648,6 +2709,7 @@ namespace
       _polygons.resize( _polygons.size() + 1 );
       _Face* polygon = &_polygons.back();
       polygon->_polyLinks.reserve( 20 );
+      polygon->_name = quad._name;
 
       splits.clear();
       for ( int iE = 0; iE < 4; ++iE ) // loop on 4 sides of a quadrangle
@@ -2682,6 +2744,7 @@ namespace
           _polygons.resize( _polygons.size() + 1 );
           polygon = &_polygons.back();
           polygon->_polyLinks.reserve( 20 );
+          polygon->_name = quad._name;
         }
         polygon->_links.push_back( splits[ iS ] );
         splits[ iS++ ]._link = 0;
@@ -2697,7 +2760,7 @@ namespace
           n1 = split.FirstNode();
           if ( n1 == n2 &&
                n1->_intPoint &&
-               (( n1->_intPoint->_faceIDs.size() > 1 && isImplementEdges() ) ||
+               (( n1->_intPoint->_faceIDs.size() > 1 && toCheckSideDivision ) ||
                 ( n1->_isInternalFlags )))
           {
             // n1 is at intersection with EDGE
@@ -3148,6 +3211,40 @@ namespace
     if ( _hasTooSmall )
       return false; // too small volume
 
+
+    // Try to find out names of no-name polygons (issue # 19887)
+    if ( _grid->IsToRemoveExcessEntities() && _polygons.back()._name == SMESH_Block::ID_NONE )
+    {
+      gp_XYZ uvwCenter =
+        0.5 * ( _grid->_coords[0][_i] + _grid->_coords[0][_i+1] ) * _grid->_axes[0] +
+        0.5 * ( _grid->_coords[1][_j] + _grid->_coords[1][_j+1] ) * _grid->_axes[1] +
+        0.5 * ( _grid->_coords[2][_k] + _grid->_coords[2][_k+1] ) * _grid->_axes[2];
+      for ( size_t i = _polygons.size() - 1; _polygons[i]._name == SMESH_Block::ID_NONE; --i )
+      {
+        _Face& face = _polygons[ i ];
+        Bnd_Box bb;
+        gp_Pnt uvw;
+        for ( size_t iL = 0; iL < face._links.size(); ++iL )
+        {
+          _Node* n = face._links[ iL ].FirstNode();
+          gp_XYZ p = SMESH_NodeXYZ( n->Node() );
+          _grid->ComputeUVW( p, uvw.ChangeCoord().ChangeData() );
+          bb.Add( uvw );
+        }
+        gp_Pnt pMin = bb.CornerMin();
+        if ( bb.IsXThin( _grid->_tol ))
+          face._name = pMin.X() < uvwCenter.X() ? SMESH_Block::ID_F0yz : SMESH_Block::ID_F1yz;
+        else if ( bb.IsYThin( _grid->_tol ))
+          face._name = pMin.Y() < uvwCenter.Y() ? SMESH_Block::ID_Fx0z : SMESH_Block::ID_Fx1z;
+        else if ( bb.IsZThin( _grid->_tol ))
+          face._name = pMin.Z() < uvwCenter.Z() ? SMESH_Block::ID_Fxy0 : SMESH_Block::ID_Fxy1;
+      }
+    }
+
+    _volumeDefs._nodes.clear();
+    _volumeDefs._quantities.clear();
+    _volumeDefs._names.clear();
+
     // create a classic cell if possible
 
     int nbPolygons = 0;
@@ -3168,14 +3265,12 @@ namespace
     else if ( nbNodes == 5 && nbPolygons == 5 ) isClassicElem = addPyra ();
     if ( !isClassicElem )
     {
-      _volumeDefs._nodes.clear();
-      _volumeDefs._quantities.clear();
-
       for ( size_t iF = 0; iF < _polygons.size(); ++iF )
       {
         const size_t nbLinks = _polygons[ iF ]._links.size();
         if ( nbLinks == 0 ) continue;
         _volumeDefs._quantities.push_back( nbLinks );
+        _volumeDefs._names.push_back( _polygons[ iF ]._name );
         for ( size_t iL = 0; iL < nbLinks; ++iL )
           _volumeDefs._nodes.push_back( _polygons[ iF ]._links[ iL ].FirstNode() );
       }
@@ -3199,7 +3294,7 @@ namespace
     int nbIntHex = 0;
 
     // set intersection nodes from GridLine's to links of allHexa
-    int i,j,k, cellIndex;
+    int i,j,k, cellIndex, iLink;
     for ( int iDir = 0; iDir < 3; ++iDir )
     {
       // loop on GridLine's parallel to iDir
@@ -3216,7 +3311,7 @@ namespace
           fourCells.Init( lineInd.I(), lineInd.J(), lineInd.K() );
           for ( int iL = 0; iL < 4; ++iL ) // loop on 4 cells sharing a link
           {
-            if ( !fourCells.GetCell( iL, i,j,k, cellIndex ))
+            if ( !fourCells.GetCell( iL, i,j,k, cellIndex, iLink ))
               continue;
             Hexahedron *& hex = allHexa[ cellIndex ];
             if ( !hex)
@@ -3224,7 +3319,6 @@ namespace
               hex = new Hexahedron( *this, i, j, k, cellIndex );
               ++nbIntHex;
             }
-            const int iLink = iL + iDir * 4;
             hex->_hexLinks[iLink]._fIntPoints.push_back( &(*ip) );
             hex->_nbFaceIntNodes += bool( ip->_node );
           }
@@ -3247,7 +3341,9 @@ namespace
       if ( hex ) // split hexahedron
       {
         intHexa.push_back( hex );
-        if ( hex->_nbFaceIntNodes > 0 || hex->_eIntPoints.size() > 0 )
+        if ( hex->_nbFaceIntNodes > 0 ||
+             hex->_eIntPoints.size() > 0 ||
+             hex->getSolids( solidIDs ) > 1 )
           continue; // treat intersected hex later in parallel
         this->init( hex->_i, hex->_j, hex->_k );
       }
@@ -3294,23 +3390,34 @@ namespace
       }
     }
 
-    // add elements resulted from hexadron intersection
+    // compute definitions of volumes resulted from hexadron intersection
 #ifdef WITH_TBB
     tbb::parallel_for ( tbb::blocked_range<size_t>( 0, intHexa.size() ),
                         ParallelHexahedron( intHexa ),
-                        tbb::simple_partitioner()); // ComputeElements() is called here
-    for ( size_t i = 0; i < intHexa.size(); ++i )
-      if ( Hexahedron * hex = intHexa[ i ] )
-        nbAdded += hex->addVolumes( helper );
+                        tbb::simple_partitioner()); // computeElements() is called here
 #else
     for ( size_t i = 0; i < intHexa.size(); ++i )
       if ( Hexahedron * hex = intHexa[ i ] )
-      {
-        hex->ComputeElements();
-        nbAdded += hex->addVolumes( helper );
-      }
+        hex->computeElements();
 #endif
 
+    // simplify polyhedrons
+    if ( _grid->IsToRemoveExcessEntities() )
+    {
+      for ( size_t i = 0; i < intHexa.size(); ++i )
+        if ( Hexahedron * hex = intHexa[ i ] )
+          hex->removeExcessSideDivision( allHexa );
+
+      for ( size_t i = 0; i < intHexa.size(); ++i )
+        if ( Hexahedron * hex = intHexa[ i ] )
+          hex->removeExcessNodes( allHexa );
+    }
+
+    // add volumes
+    for ( size_t i = 0; i < intHexa.size(); ++i )
+      if ( Hexahedron * hex = intHexa[ i ] )
+        nbAdded += hex->addVolumes( helper );
+
     // fill boundaryVolumes with volumes neighboring too small skipped volumes
     if ( _grid->_toCreateFaces )
     {
@@ -3604,7 +3711,8 @@ namespace
           int  i = ! ( u < _grid->_tol ); // [0,1]
           int iN = link._nodes[ i ] - hex->_hexNodes; // [0-7]
 
-          const F_IntersectPoint * & ip = _grid->_gridIntP[ hex->_origNodeInd + _nodeShift[iN] ];
+          const F_IntersectPoint * & ip = _grid->_gridIntP[ hex->_origNodeInd +
+                                                            _grid->_nodeShift[iN] ];
           if ( !ip )
           {
             ip = _grid->_extIntPool.getNew();
@@ -3632,13 +3740,12 @@ namespace
           int i,j,k, cellIndex;
           for ( int iC = 0; iC < 4; ++iC ) // loop on 4 cells sharing the link
           {
-            if ( !fourCells.GetCell( iC, i,j,k, cellIndex ))
+            if ( !fourCells.GetCell( iC, i,j,k, cellIndex, iLink ))
               continue;
             Hexahedron * h = hexes[ cellIndex ];
             if ( !h )
               h = hexes[ cellIndex ] = new Hexahedron( *this, i, j, k, cellIndex );
-            const int iL = iC + iDir * 4;
-            h->_hexLinks[iL]._fIntPoints.push_back( ip );
+            h->_hexLinks[iLink]._fIntPoints.push_back( ip );
             h->_nbFaceIntNodes++;
             //isCut = true;
           }
@@ -3698,8 +3805,8 @@ namespace
 
     for ( size_t iN = 0; iN < 8; ++iN ) // check corners
     {
-      _hexNodes[iN]._node     = _grid->_nodes   [ _origNodeInd + _nodeShift[iN] ];
-      _hexNodes[iN]._intPoint = _grid->_gridIntP[ _origNodeInd + _nodeShift[iN] ];
+      _hexNodes[iN]._node     = _grid->_nodes   [ _origNodeInd + _grid->_nodeShift[iN] ];
+      _hexNodes[iN]._intPoint = _grid->_gridIntP[ _origNodeInd + _grid->_nodeShift[iN] ];
       if ( _hexNodes[iN]._intPoint )
         for ( size_t iF = 0; iF < _hexNodes[iN]._intPoint->_faceIDs.size(); ++iF )
         {
@@ -4369,10 +4476,13 @@ namespace
         {
           curIntPnt._paramOnLine = coords[ ijk[ iDir ]] - coords[0] + _grid->_tol;
           const GridLine& line = _grid->_lines[ iDir ][ lineIndex[ iL ]];
-          multiset< F_IntersectPoint >::const_iterator ip =
-            line._intPoints.upper_bound( curIntPnt );
-          --ip;
-          firstIntPnt = &(*ip);
+          if ( !line._intPoints.empty() )
+          {
+            multiset< F_IntersectPoint >::const_iterator ip =
+              line._intPoints.upper_bound( curIntPnt );
+            --ip;
+            firstIntPnt = &(*ip);
+          }
         }
         else if ( !link._fIntPoints.empty() )
         {
@@ -4958,7 +5068,7 @@ namespace
   /*!
    * \brief Return created volumes and volumes that can have free facet because of
    *        skipped small volume. Also create mesh faces on free facets
-   *        of adjacent not-cut volumes id the result volume is too small.
+   *        of adjacent not-cut volumes if the result volume is too small.
    */
   void Hexahedron::getBoundaryElems( vector< const SMDS_MeshElement* > & boundaryElems )
   {
@@ -5005,8 +5115,9 @@ namespace
           if ( !faceID )
             break;
           if ( _grid->IsInternal( faceID ) ||
-               _grid->IsShared( faceID ) ||
-               _grid->IsBoundaryFace( faceID ))
+               _grid->IsShared( faceID ) //||
+               //_grid->IsBoundaryFace( faceID ) -- commented for #19887
+               ) 
             break; // create only if a new face will be used by other 3D algo
         }
 
@@ -5035,6 +5146,400 @@ namespace
     }
   }
 
+  //================================================================================
+  /*!
+   * \brief Remove edges and nodes dividing a hexa side in the case if an adjacent
+   *        volume also sharing the dividing edge is missing due to its small side.
+   *        Issue #19887.
+   */
+  //================================================================================
+
+  void Hexahedron::removeExcessSideDivision(const vector< Hexahedron* >& allHexa)
+  {
+    if ( ! _volumeDefs.IsPolyhedron() )
+      return; // not a polyhedron
+      
+    // look for a divided side adjacent to a small hexahedron
+
+    int di[6] = { 0, 0, 0, 0,-1, 1 };
+    int dj[6] = { 0, 0,-1, 1, 0, 0 };
+    int dk[6] = {-1, 1, 0, 0, 0, 0 };
+
+    for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides of a hexahedron
+    {
+      size_t neighborIndex = _grid->CellIndex( _i + di[iF],
+                                               _j + dj[iF],
+                                               _k + dk[iF] );
+      if ( neighborIndex >= allHexa.size() ||
+           !allHexa[ neighborIndex ]       ||
+           !allHexa[ neighborIndex ]->_hasTooSmall )
+        continue;
+
+      // check if a side is divided into several polygons
+      for ( _volumeDef* volDef = &_volumeDefs; volDef; volDef = volDef->_next )
+      {
+        int nbPolygons = 0, nbNodes = 0;
+        for ( size_t i = 0; i < volDef->_names.size(); ++i )
+          if ( volDef->_names[ i ] == _hexQuads[ iF ]._name )
+          {
+            ++nbPolygons;
+            nbNodes += volDef->_quantities[ i ];
+          }
+        if ( nbPolygons < 2 )
+          continue;
+
+        // construct loops from polygons
+        typedef _volumeDef::_linkDef TLinkDef;
+        std::vector< TLinkDef* > loops;
+        std::vector< TLinkDef > links( nbNodes );
+        for ( size_t i = 0, iN = 0, iLoop = 0; iLoop < volDef->_quantities.size(); ++iLoop )
+        {
+          size_t nbLinks = volDef->_quantities[ iLoop ];
+          if ( volDef->_names[ iLoop ] != _hexQuads[ iF ]._name )
+          {
+            iN += nbLinks;
+            continue;
+          }
+          loops.push_back( & links[i] );
+          for ( size_t n = 0; n < nbLinks-1; ++n, ++i, ++iN )
+          {
+            links[i].init( volDef->_nodes[iN], volDef->_nodes[iN+1], iLoop );
+            links[i].setNext( &links[i+1] );
+          }
+          links[i].init( volDef->_nodes[iN], volDef->_nodes[iN-nbLinks+1], iLoop );
+          links[i].setNext( &links[i-nbLinks+1] );
+          ++i; ++iN;
+        }
+
+        // look for equal links in different loops and join such loops
+        bool loopsJoined = false;
+        std::set< TLinkDef > linkSet;
+        for ( size_t iLoop = 0; iLoop < loops.size(); ++iLoop )
+        {
+          TLinkDef* beg = 0;
+          for ( TLinkDef* l = loops[ iLoop ]; l != beg; l = l->_next ) // walk around the iLoop
+          {
+            std::pair< std::set< TLinkDef >::iterator, bool > it2new = linkSet.insert( *l );
+            if ( !it2new.second ) // equal found, join loops
+            {
+              const TLinkDef* equal = &(*it2new.first);
+              if ( equal->_loopIndex == l->_loopIndex )
+                continue; // error?
+
+              loopsJoined = true;
+
+              for ( size_t i = iLoop - 1; i < loops.size(); --i )
+                if ( loops[ i ] && loops[ i ]->_loopIndex == equal->_loopIndex )
+                  loops[ i ] = 0;
+
+              // exclude l and equal and join two loops
+              if ( l->_prev != equal )
+                l->_prev->setNext( equal->_next );
+              if ( equal->_prev != l )
+                equal->_prev->setNext( l->_next );
+
+              if ( volDef->_quantities[ l->_loopIndex ] > 0 )
+                volDef->_quantities[ l->_loopIndex     ] *= -1;
+              if ( volDef->_quantities[ equal->_loopIndex ] > 0 )
+                volDef->_quantities[ equal->_loopIndex ] *= -1;
+
+              if ( loops[ iLoop ] == l )
+                loops[ iLoop ] = l->_prev->_next;
+            }
+            beg = loops[ iLoop ];
+          }
+        }
+        // update volDef
+        if ( loopsJoined )
+        {
+          // set unchanged polygons
+          std::vector< int >                  newQuantities;
+          std::vector< _volumeDef::_nodeDef > newNodes;
+          vector< SMESH_Block::TShapeID >     newNames;
+          newQuantities.reserve( volDef->_quantities.size() );
+          newNodes.reserve     ( volDef->_nodes.size() );
+          newNames.reserve     ( volDef->_names.size() );
+          for ( size_t i = 0, iLoop = 0; iLoop < volDef->_quantities.size(); ++iLoop )
+          {
+            if ( volDef->_quantities[ iLoop ] < 0 )
+            {
+              i -= volDef->_quantities[ iLoop ];
+              continue;
+            }
+            newQuantities.push_back( volDef->_quantities[ iLoop ]);
+            newNodes.insert( newNodes.end(),
+                             volDef->_nodes.begin() + i,
+                             volDef->_nodes.begin() + i + newQuantities.back() );
+            newNames.push_back( volDef->_names[ iLoop ]);
+            i += volDef->_quantities[ iLoop ];
+          }
+
+          // set joined loops
+          for ( size_t iLoop = 0; iLoop < loops.size(); ++iLoop )
+          {
+            if ( !loops[ iLoop ] )
+              continue;
+            newQuantities.push_back( 0 );
+            TLinkDef* beg = 0;
+            for ( TLinkDef* l = loops[ iLoop ]; l != beg; l = l->_next, ++newQuantities.back() )
+            {
+              newNodes.push_back( l->_node1 );
+              beg = loops[ iLoop ];
+            }
+            newNames.push_back( _hexQuads[ iF ]._name );
+          }
+          volDef->_quantities.swap( newQuantities );
+          volDef->_nodes.swap( newNodes );
+          volDef->_names.swap( newNames );
+        }
+      } // loop on volDef's
+    } // loop on hex sides
+
+    return;
+  } // removeExcessSideDivision()
+
+
+  //================================================================================
+  /*!
+   * \brief Remove nodes splitting Cartesian cell edges in the case if a node
+   *        is used in every cells only by two polygons sharing the edge
+   *        Issue #19887.
+   */
+  //================================================================================
+
+  void Hexahedron::removeExcessNodes(vector< Hexahedron* >& allHexa)
+  {
+    if ( ! _volumeDefs.IsPolyhedron() )
+      return; // not a polyhedron
+
+    typedef vector< _volumeDef::_nodeDef >::iterator TNodeIt;
+    vector< int > nodesInPoly[ 4 ]; // node index in _volumeDefs._nodes
+    vector< int > volDefInd  [ 4 ]; // index of a _volumeDefs
+    Hexahedron*   hexa       [ 4 ];
+    int i,j,k, cellIndex, iLink = 0, iCellLink;
+    for ( int iDir = 0; iDir < 3; ++iDir )
+    {
+      CellsAroundLink fourCells( _grid, iDir );
+      for ( int iL = 0; iL < 4; ++iL, ++iLink ) // 4 links in a direction
+      {
+        _Link& link = _hexLinks[ iLink ];
+        fourCells.Init( _i, _j, _k, iLink );
+
+        for ( size_t iP = 0; iP < link._fIntPoints.size(); ++iP ) // loop on nodes on the link
+        {
+          bool nodeRemoved = true;
+          _volumeDef::_nodeDef node; node._intPoint = link._fIntPoints[iP];
+
+          for ( size_t i = 0, nb = _volumeDefs.size(); i < nb &&  nodeRemoved; ++i )
+            if ( _volumeDef* vol = _volumeDefs.at( i ))
+              nodeRemoved =
+                ( std::find( vol->_nodes.begin(), vol->_nodes.end(), node ) == vol->_nodes.end() );
+          if ( nodeRemoved )
+            continue; // node already removed
+
+          // check if a node encounters zero or two times in 4 cells sharing iLink
+          // if so, the node can be removed from the cells
+          bool       nodeIsOnEdge = true;
+          int nbPolyhedraWithNode = 0;
+          for ( int iC = 0; iC < 4; ++iC ) // loop on 4 cells sharing a link
+          {
+            nodesInPoly[ iC ].clear();
+            volDefInd  [ iC ].clear();
+            hexa       [ iC ] = 0;
+            if ( !fourCells.GetCell( iC, i,j,k, cellIndex, iCellLink ))
+              continue;
+            hexa[ iC ] = allHexa[ cellIndex ];
+            if ( !hexa[ iC ])
+              continue;
+            for ( size_t i = 0, nb = hexa[ iC ]->_volumeDefs.size(); i < nb; ++i )
+              if ( _volumeDef* vol = hexa[ iC ]->_volumeDefs.at( i ))
+              {
+                for ( TNodeIt nIt = vol->_nodes.begin(); nIt != vol->_nodes.end(); ++nIt )
+                {
+                  nIt = std::find( nIt, vol->_nodes.end(), node );
+                  if ( nIt != vol->_nodes.end() )
+                  {
+                    nodesInPoly[ iC ].push_back( std::distance( vol->_nodes.begin(), nIt ));
+                    volDefInd  [ iC ].push_back( i );
+                  }
+                  else
+                    break;
+                }
+                nbPolyhedraWithNode += ( !nodesInPoly[ iC ].empty() );
+              }
+            if ( nodesInPoly[ iC ].size() != 0 &&
+                 nodesInPoly[ iC ].size() != 2 )
+            {
+              nodeIsOnEdge = false;
+              break;
+            }
+          } // loop  on 4 cells
+
+          // remove nodes from polyhedra
+          if ( nbPolyhedraWithNode > 0 && nodeIsOnEdge )
+          {
+            for ( int iC = 0; iC < 4; ++iC ) // loop on 4 cells sharing the link
+            {
+              if ( nodesInPoly[ iC ].empty() )
+                continue;
+              for ( int i = volDefInd[ iC ].size() - 1; i >= 0; --i )
+              {
+                _volumeDef* vol = hexa[ iC ]->_volumeDefs.at( volDefInd[ iC ][ i ]);
+                int nIndex = nodesInPoly[ iC ][ i ];
+                // decrement _quantities
+                for ( size_t iQ = 0; iQ < vol->_quantities.size(); ++iQ )
+                  if ( nIndex < vol->_quantities[ iQ ])
+                  {
+                    vol->_quantities[ iQ ]--;
+                    break;
+                  }
+                  else
+                  {
+                    nIndex -= vol->_quantities[ iQ ];
+                  }
+                vol->_nodes.erase( vol->_nodes.begin() + nodesInPoly[ iC ][ i ]);
+
+                if ( i == 0 &&
+                     vol->_nodes.size() == 6 * 4 &&
+                     vol->_quantities.size() == 6 ) // polyhedron becomes hexahedron?
+                {
+                  bool allQuads = true;
+                  for ( size_t iQ = 0; iQ < vol->_quantities.size() &&  allQuads; ++iQ )
+                    allQuads = ( vol->_quantities[ iQ ] == 4 );
+                  if ( allQuads )
+                  {
+                    // set side nodes as this: bottom, top, top, ...
+                    int iTop, iBot; // side indices
+                    for ( int iS = 0; iS < 6; ++iS )
+                    {
+                      if ( vol->_names[ iS ] == SMESH_Block::ID_Fxy0 )
+                        iBot = iS;
+                      if ( vol->_names[ iS ] == SMESH_Block::ID_Fxy1 )
+                        iTop = iS;
+                    }
+                    if ( iBot != 0 )
+                    {
+                      if ( iTop == 0 )
+                      {
+                        std::copy( vol->_nodes.begin(),
+                                   vol->_nodes.begin() + 4,
+                                   vol->_nodes.begin() + 4 );
+                        iTop = 1;
+                      }
+                      std::copy( vol->_nodes.begin() + 4 * iBot,
+                                 vol->_nodes.begin() + 4 * ( iBot + 1),
+                                 vol->_nodes.begin() );
+                    }
+                    if ( iTop != 1 )
+                      std::copy( vol->_nodes.begin() + 4 * iTop,
+                                 vol->_nodes.begin() + 4 * ( iTop + 1),
+                                 vol->_nodes.begin() + 4 );
+
+                    std::copy( vol->_nodes.begin() + 4,
+                               vol->_nodes.begin() + 8,
+                               vol->_nodes.begin() + 8 );
+                    // set up top facet nodes by comparing their uvw with bottom nodes
+                    E_IntersectPoint ip[8];
+                    for ( int iN = 0; iN < 8; ++iN )
+                    {
+                      SMESH_NodeXYZ p = vol->_nodes[ iN ].Node();
+                      _grid->ComputeUVW( p, ip[ iN ]._uvw );
+                    }
+                    const double tol2 = _grid->_tol * _grid->_tol;
+                    for ( int iN = 0; iN < 4; ++iN )
+                    {
+                      gp_Pnt2d pBot( ip[ iN ]._uvw[0], ip[ iN ]._uvw[1] );
+                      for ( int iT = 4; iT < 8; ++iT )
+                      {
+                        gp_Pnt2d pTop( ip[ iT ]._uvw[0], ip[ iT ]._uvw[1] );
+                        if ( pBot.SquareDistance( pTop ) < tol2 )
+                        {
+                          // vol->_nodes[ iN + 4 ]._node = ip[ iT ]._node;
+                          // vol->_nodes[ iN + 4 ]._intPoint = 0;
+                          vol->_nodes[ iN + 4 ] = vol->_nodes[ iT + 4 ];
+                          break;
+                        }
+                      }
+                    }
+                    vol->_nodes.resize( 8 );
+                    vol->_quantities.clear();
+                    //vol->_names.clear();
+                  }
+                }
+              } // loop on _volumeDefs
+            } // loop on 4 cell abound a link
+          } // if ( nodeIsOnEdge )
+        } // loop on intersection points of a link
+      } // loop on 4 links of a direction
+    } // loop on 3 directions
+
+    return;
+
+  } // removeExcessNodes()
+
+  //================================================================================
+  /*!
+   * \brief [Issue #19913] Modify _hexLinks._splits to prevent creating overlapping volumes
+   */
+  //================================================================================
+
+  void Hexahedron::preventVolumesOverlapping()
+  {
+    // Cut off a quadrangle corner if two links sharing the corner
+    // are shared by same two solids, in this case each of solids gets
+    // a triangle for it-self.
+    std::vector< TGeomID > soIDs[4];
+    for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides of a hexahedron
+    {
+      _Face& quad = _hexQuads[ iF ] ;
+
+      int iFOpposite = iF + ( iF % 2 ? -1 : 1 );
+      _Face& quadOpp = _hexQuads[ iFOpposite ] ;
+
+      int nbSides = 0, nbSidesOpp = 0;
+      for ( int iE = 0; iE < 4; ++iE ) // loop on 4 sides of a quadrangle
+      {
+        nbSides    += ( quad._links   [ iE ].NbResultLinks() > 0 );
+        nbSidesOpp += ( quadOpp._links[ iE ].NbResultLinks() > 0 );
+      }
+      if ( nbSides < 4 || nbSidesOpp != 2 )
+        continue;
+
+      for ( int iE = 0; iE < 4; ++iE )
+      {
+        soIDs[ iE ].clear();
+        _Node* n = quad._links[ iE ].FirstNode();
+        if ( n->_intPoint && n->_intPoint->_faceIDs.size() )
+          soIDs[ iE ] = _grid->GetSolidIDs( n->_intPoint->_faceIDs[0] );
+      }
+      if ((( soIDs[0].size() >= 2 ) +
+           ( soIDs[1].size() >= 2 ) +
+           ( soIDs[2].size() >= 2 ) +
+           ( soIDs[3].size() >= 2 ) ) < 3 )
+        continue;
+
+      bool done = false;
+      for ( int i = 0; i < 4; ++i )
+      {
+        int i1 = _grid->_helper->WrapIndex( i + 1, 4 );
+        int i2 = _grid->_helper->WrapIndex( i + 2, 4 );
+        int i3 = _grid->_helper->WrapIndex( i + 3, 4 );
+        if ( soIDs[i1].size() == 2 && soIDs[i ] != soIDs[i1] &&
+             soIDs[i2].size() == 2 && soIDs[i1] == soIDs[i2] &&
+             soIDs[i3].size() == 2 && soIDs[i2] == soIDs[i3] )
+        {
+          quad._links[ i1 ]._link->_splits.clear();
+          quad._links[ i2 ]._link->_splits.clear();
+          done = true;
+          break;
+        }
+      }
+      if ( done )
+        break;
+    }
+    return;
+  } // preventVolumesOverlapping()
+
   //================================================================================
   /*!
    * \brief Set to _hexLinks a next portion of splits located on one side of INTERNAL FACEs
diff --git a/src/StdMeshers/StdMeshers_Cartesian_3D.hxx b/src/StdMeshers/StdMeshers_Cartesian_3D.hxx
index 593646a74..18797e278 100644
--- a/src/StdMeshers/StdMeshers_Cartesian_3D.hxx
+++ b/src/StdMeshers/StdMeshers_Cartesian_3D.hxx
@@ -35,6 +35,8 @@
  * internal part of a solid shape and polyhedral volumes near the shape boundary.
  *
  * Issue 0021336
+ * Issue #16523: Treatment of internal faces
+ * Issue #17237: Body fitting on sub-mesh
  */
 class StdMeshers_CartesianParameters3D;
 
diff --git a/src/StdMeshersGUI/StdMeshersGUI_SubShapeSelectorWdg.cxx b/src/StdMeshersGUI/StdMeshersGUI_SubShapeSelectorWdg.cxx
index 54560d8fc..31db0a7be 100644
--- a/src/StdMeshersGUI/StdMeshersGUI_SubShapeSelectorWdg.cxx
+++ b/src/StdMeshersGUI/StdMeshersGUI_SubShapeSelectorWdg.cxx
@@ -74,7 +74,8 @@ StdMeshersGUI_SubShapeSelectorWdg
 ::StdMeshersGUI_SubShapeSelectorWdg( QWidget *        parent,
                                      TopAbs_ShapeEnum subShType,
                                      const bool       toShowList,
-                                     const bool       toShowActivateBtn ):
+                                     const bool       toShowActivateBtn,
+                                     const int        minListWidth ):
   QWidget( parent ),
   myMaxSize( -1 ),
   myPreviewActor( 0 )
@@ -93,8 +94,8 @@ StdMeshersGUI_SubShapeSelectorWdg
     myAddButton      = new QPushButton( tr( "SMESH_BUT_ADD" ),    this );
     myRemoveButton   = new QPushButton( tr( "SMESH_BUT_REMOVE" ), this );
     myListWidget->setSelectionMode( QListWidget::ExtendedSelection );
-    myListWidget->setMinimumWidth(300);
-    myListWidget->setWrapping(true);
+    myListWidget->setMinimumWidth( minListWidth );
+    myListWidget->setWrapping( true );
     myActivateButton->setCheckable( true );
   }
   else
@@ -254,6 +255,20 @@ void StdMeshersGUI_SubShapeSelectorWdg::ShowPreview( bool visible)
   }
 }
 
+//================================================================================
+/*!
+ * \brief Connect selection slots
+  * \param other - another StdMeshersGUI_ObjectReferenceParamWdg
+ */
+//================================================================================
+
+void StdMeshersGUI_SubShapeSelectorWdg::
+AvoidSimultaneousSelection ( StdMeshersGUI_SubShapeSelectorWdg* other)
+{
+  connect(other, SIGNAL(selectionActivated()), this, SLOT(deactivateSelection()));
+  connect(this, SIGNAL(selectionActivated()), other, SLOT(deactivateSelection()));
+}
+
 //================================================================================
 /*!
  * \brief Connect/disconnect to change of selection
@@ -275,7 +290,7 @@ void StdMeshersGUI_SubShapeSelectorWdg::ActivateSelection( bool toActivate )
 
   if ( toActivate )
   {
-    connect( mySelectionMgr, SIGNAL(currentSelectionChanged()), this, SLOT(selectionIntoArgument()));
+    connect( mySelectionMgr, SIGNAL(currentSelectionChanged()), SLOT(selectionIntoArgument()));
   }
   else
   {
@@ -284,6 +299,9 @@ void StdMeshersGUI_SubShapeSelectorWdg::ActivateSelection( bool toActivate )
 
   if ( sender() == myActivateButton )
     ShowPreview( toActivate );
+
+  if ( toActivate )
+    emit selectionActivated();
 }
 
 //================================================================================
diff --git a/src/StdMeshersGUI/StdMeshersGUI_SubShapeSelectorWdg.h b/src/StdMeshersGUI/StdMeshersGUI_SubShapeSelectorWdg.h
index 3f3d3bf12..9c4f14b35 100644
--- a/src/StdMeshersGUI/StdMeshersGUI_SubShapeSelectorWdg.h
+++ b/src/StdMeshersGUI/StdMeshersGUI_SubShapeSelectorWdg.h
@@ -53,7 +53,8 @@ public:
   StdMeshersGUI_SubShapeSelectorWdg( QWidget*         parent = 0,
                                      TopAbs_ShapeEnum subShType = TopAbs_EDGE,
                                      const bool       toShowList = true,
-                                     const bool       toShowSelectBtn = false);
+                                     const bool       toShowSelectBtn = false,
+                                     const int        minListWidth=300);
   ~StdMeshersGUI_SubShapeSelectorWdg();
 
   SMESH::long_array_var          GetListOfIDs();
@@ -82,12 +83,23 @@ public:
   SMESH_PreviewActorsCollection* GetActorCollection() { return myPreviewActor; }
   void                           ClearSelected();
 
+  void AvoidSimultaneousSelection( StdMeshersGUI_SubShapeSelectorWdg* other);
+
 public slots:
+
   void                           ActivateSelection( bool );
+  void                           deactivateSelection() { ActivateSelection( false ); }
 
 signals:
   void                           selectionChanged(); // in the list
   void                           shapeSelected();    // globally
+  /*!
+   * \brief Emitted when selection is activated
+    *
+    * Useful to deactivate one Object Reference param widget when an other
+    * one is activated
+   */
+  void                           selectionActivated();
 
 private:
   void                           updateState();