Salome HOME
remove trace
[modules/geom.git] / src / GEOMImpl / GEOMImpl_IBlocksOperations.cxx
index 9b7aedc80d8bb867edd7ef7d32697ddf0c56c687..b08d2b829e9ca7b2a275092da4ffd3f6fb3900f5 100644 (file)
@@ -1,28 +1,67 @@
-using namespace std;
-
-#include "GEOMImpl_IBlocksOperations.hxx"
-
-#include "GEOMImpl_Types.hxx"
-
-#include "GEOMImpl_BlockDriver.hxx"
-#include "GEOMImpl_IBlocks.hxx"
-#include "GEOMImpl_IBlockTrsf.hxx"
-#include "GEOMImpl_CopyDriver.hxx"
-#include "GEOMImpl_Block6Explorer.hxx"
-
-#include "GEOM_Function.hxx"
+// Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifdef WIN32
+#pragma warning( disable:4786 )
+#endif
+
+#include <Standard_Stream.hxx>
+
+#include <GEOMImpl_IBlocksOperations.hxx>
+
+#include <GEOMImpl_Types.hxx>
+
+#include <GEOMImpl_BlockDriver.hxx>
+#include <GEOMImpl_IBlocks.hxx>
+#include <GEOMImpl_IBlockTrsf.hxx>
+#include <GEOMImpl_Block6Explorer.hxx>
+
+#include <GEOMUtils.hxx>
+
+#include <GEOM_Function.hxx>
+#include <GEOM_PythonDump.hxx>
+
+#include <GEOMAlgo_GlueAnalyser.hxx>
+#include <GEOMAlgo_CoupleOfShapes.hxx>
+#include <GEOMAlgo_ListOfCoupleOfShapes.hxx>
+#include <GEOMAlgo_ListIteratorOfListOfCoupleOfShapes.hxx>
+#include <BlockFix_CheckTool.hxx>
+
+#include <Basics_OCCTVersion.hxx>
 
 #include "utilities.h"
-#include "OpUtil.hxx"
-#include "Utils_ExceptHandlers.hxx"
+#include <OpUtil.hxx>
+#include <Utils_ExceptHandlers.hxx>
 
 #include <TFunction_DriverTable.hxx>
 #include <TFunction_Driver.hxx>
 #include <TFunction_Logbook.hxx>
+#include <TDataStd_Integer.hxx>
 #include <TDF_Tool.hxx>
 
 #include <BRep_Tool.hxx>
+#include <BRep_Builder.hxx>
 #include <BRepTools.hxx>
+#include <BRepTools_WireExplorer.hxx>
 #include <BRepGProp.hxx>
 #include <BRepBndLib.hxx>
 #include <BRepAdaptor_Surface.hxx>
@@ -34,23 +73,34 @@ using namespace std;
 #include <TopoDS.hxx>
 #include <TopoDS_Edge.hxx>
 #include <TopoDS_Vertex.hxx>
+#include <TopoDS_Compound.hxx>
 #include <TopoDS_Iterator.hxx>
 #include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
 #include <TopTools_MapOfShape.hxx>
 #include <TopTools_Array1OfShape.hxx>
 #include <TopTools_IndexedMapOfShape.hxx>
+#include <TopTools_DataMapOfShapeInteger.hxx>
 #include <TopTools_ListIteratorOfListOfShape.hxx>
 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
+#include <TopTools_DataMapIteratorOfDataMapOfShapeInteger.hxx>
 
 #include <Bnd_Box.hxx>
-#include <Precision.hxx>
 #include <GProp_GProps.hxx>
+
+#include <Geom_Surface.hxx>
+#include <ShapeAnalysis_Surface.hxx>
+
 #include <TColStd_MapOfInteger.hxx>
 #include <TColStd_Array1OfReal.hxx>
 #include <TColStd_Array1OfInteger.hxx>
 #include <TColStd_Array2OfInteger.hxx>
 
+//#include <OSD_Timer.hxx>
+
+#include <Precision.hxx>
+
+#include <Standard_Failure.hxx>
 #include <Standard_ErrorHandler.hxx> // CAREFUL ! position of this file is critic : see Lucien PIGNOLONI / OCC
 
 //=============================================================================
@@ -118,6 +168,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad
 
   //Compute the Face value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to compute a face");
       return NULL;
@@ -130,19 +183,8 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad
   }
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aFace->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.MakeQuad(";
-  TDF_Tool::Entry(theEdge1->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(theEdge2->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(theEdge3->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(theEdge4->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  aFunction->SetDescription(aDescr);
+  GEOM::TPythonDump(aFunction) << aFace << " = geompy.MakeQuad("
+    << theEdge1 << ", " << theEdge2 << ", " << theEdge3 << ", " << theEdge4 << ")";
 
   SetErrorCode(OK);
   return aFace;
@@ -184,6 +226,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad2Edges
 
   //Compute the Face value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to compute a face");
       return NULL;
@@ -196,15 +241,8 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad2Edges
   }
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aFace->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.MakeQuad2Edges(";
-  TDF_Tool::Entry(theEdge1->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(theEdge2->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  aFunction->SetDescription(aDescr);
+  GEOM::TPythonDump(aFunction) << aFace << " = geompy.MakeQuad2Edges("
+                               << theEdge1 << ", " << theEdge2 << ")";
 
   SetErrorCode(OK);
   return aFace;
@@ -253,6 +291,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad4Vertices
 
   //Compute the Face value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to compute a face");
       return NULL;
@@ -265,19 +306,8 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad4Vertices
   }
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aFace->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.MakeQuad4Vertices(";
-  TDF_Tool::Entry(thePnt1->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(thePnt2->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(thePnt3->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(thePnt4->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  aFunction->SetDescription(aDescr);
+  GEOM::TPythonDump(aFunction) << aFace << " = geompy.MakeQuad4Vertices("
+    << thePnt1 << ", " << thePnt2 << ", " << thePnt3 << ", " << thePnt4 << ")";
 
   SetErrorCode(OK);
   return aFace;
@@ -333,6 +363,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeHexa
 
   //Compute the Block value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to compute a block");
       return NULL;
@@ -345,23 +378,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeHexa
   }
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aBlock->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.MakeHexa(";
-  TDF_Tool::Entry(theFace1->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(theFace2->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(theFace3->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(theFace4->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(theFace5->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(theFace6->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  aFunction->SetDescription(aDescr);
+  GEOM::TPythonDump(aFunction) << aBlock << " = geompy.MakeHexa("
+    << theFace1 << ", " << theFace2 << ", " << theFace3 << ", "
+      << theFace4 << ", " << theFace5 << ", " << theFace6 << ")";
 
   SetErrorCode(OK);
   return aBlock;
@@ -403,6 +422,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeHexa2Faces
 
   //Compute the Block value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to compute a block");
       return NULL;
@@ -415,15 +437,8 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeHexa2Faces
   }
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aBlock->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.MakeHexa2Faces(";
-  TDF_Tool::Entry(theFace1->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(theFace2->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  aFunction->SetDescription(aDescr);
+  GEOM::TPythonDump(aFunction) << aBlock << " = geompy.MakeHexa2Faces("
+                               << theFace1 << ", " << theFace2 << ")";
 
   SetErrorCode(OK);
   return aBlock;
@@ -463,6 +478,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeBlockCompound
 
   //Compute the Blocks Compound value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to compute a blocks compound");
       return NULL;
@@ -475,13 +493,8 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeBlockCompound
   }
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aBlockComp->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.MakeBlockCompound(";
-  TDF_Tool::Entry(theCompound->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  aFunction->SetDescription(aDescr);
+  GEOM::TPythonDump(aFunction) << aBlockComp
+    << " = geompy.MakeBlockCompound(" << theCompound << ")";
 
   SetErrorCode(OK);
   return aBlockComp;
@@ -489,7 +502,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeBlockCompound
 
 //=============================================================================
 /*!
- *  GetEdge
+ *  GetPoint
  */
 //=============================================================================
 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetPoint
@@ -509,13 +522,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetPoint
 
   TopoDS_Shape aBlockOrComp = theShape->GetValue();
   if (aBlockOrComp.IsNull()) {
-    SetErrorCode("Block or compound is null");
-    return NULL;
-  }
-  if (aBlockOrComp.ShapeType() != TopAbs_SOLID &&
-      aBlockOrComp.ShapeType() != TopAbs_COMPOUND &&
-      aBlockOrComp.ShapeType() != TopAbs_COMPSOLID) {
-    SetErrorCode("Shape is neither a block, nor a compound of blocks");
+    SetErrorCode("Given shape is null");
     return NULL;
   }
 
@@ -557,19 +564,82 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetPoint
   Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aResult->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.GetPoint(";
-  TDF_Tool::Entry(theShape->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  aDescr += TCollection_AsciiString(theX) + ", ";
-  aDescr += TCollection_AsciiString(theY) + ", ";
-  aDescr += TCollection_AsciiString(theZ) + ", ";
-  aDescr += TCollection_AsciiString(theEpsilon) + ")";
-
-  TCollection_AsciiString aNewDescr = aFunction->GetDescription() + "\n";
-  aNewDescr += aDescr;
-  aFunction->SetDescription(aNewDescr);
+  GEOM::TPythonDump(aFunction, /*append=*/true)
+    << aResult << " = geompy.GetPoint(" << theShape << ", "
+    << theX << ", " << theY << ", " << theZ << ", " << theEpsilon << ")";
+
+  SetErrorCode(OK);
+  return aResult;
+}
+
+//=============================================================================
+/*!
+ *  GetVertexNearPoint
+ */
+//=============================================================================
+Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetVertexNearPoint
+                                               (Handle(GEOM_Object) theShape,
+                                                Handle(GEOM_Object) thePoint)
+{
+  SetErrorCode(KO);
+
+  // New Point object
+  Handle(GEOM_Object) aResult;
+
+  // Arguments
+  if (theShape.IsNull() || thePoint.IsNull()) return NULL;
+
+  TopoDS_Shape aBlockOrComp = theShape->GetValue();
+  TopoDS_Shape aPoint       = thePoint->GetValue();
+  if (aBlockOrComp.IsNull() || aPoint.IsNull()) {
+    SetErrorCode("Given shape is null");
+    return NULL;
+  }
+
+  if (aPoint.ShapeType() != TopAbs_VERTEX) {
+    SetErrorCode("Element for vertex identification is not a vertex");
+    return NULL;
+  }
+
+  TopoDS_Vertex aVert = TopoDS::Vertex(aPoint);
+  gp_Pnt aP = BRep_Tool::Pnt(aVert);
+
+  // Compute the Vertex value
+  TopoDS_Shape V;
+  bool isFound = false;
+  Standard_Real aDist = RealLast();
+  TopTools_MapOfShape mapShape;
+
+  TopExp_Explorer exp (aBlockOrComp, TopAbs_VERTEX);
+  for (; exp.More(); exp.Next()) {
+    if (mapShape.Add(exp.Current())) {
+      TopoDS_Vertex aVi = TopoDS::Vertex(exp.Current());
+      gp_Pnt aPi = BRep_Tool::Pnt(aVi);
+      Standard_Real aDisti = aPi.Distance(aP);
+      if (aDisti < aDist) {
+        V = aVi;
+        aDist = aDisti;
+        isFound = true;
+      }
+    }
+  }
+
+  if (!isFound) {
+    SetErrorCode("Vertex has not been found");
+    return NULL;
+  }
+
+  TopTools_IndexedMapOfShape anIndices;
+  TopExp::MapShapes(aBlockOrComp, anIndices);
+  Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
+  anArray->SetValue(1, anIndices.FindIndex(V));
+  aResult = GetEngine()->AddSubShape(theShape, anArray);
+
+  Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
+
+  // Make a Python command
+  GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetVertexNearPoint("
+                               << theShape << ", " << thePoint << ")";
 
   SetErrorCode(OK);
   return aResult;
@@ -595,13 +665,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetEdge
 
   TopoDS_Shape aBlockOrComp = theShape->GetValue();
   if (aBlockOrComp.IsNull()) {
-    SetErrorCode("Block or compound is null");
-    return NULL;
-  }
-  if (aBlockOrComp.ShapeType() != TopAbs_SOLID &&
-      aBlockOrComp.ShapeType() != TopAbs_COMPOUND &&
-      aBlockOrComp.ShapeType() != TopAbs_COMPSOLID) {
-    SetErrorCode("Shape is neither a block, nor a compound of blocks");
+    SetErrorCode("Given shape is null");
     return NULL;
   }
 
@@ -619,6 +683,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetEdge
 
   //Compute the Edge value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     TopTools_IndexedDataMapOfShapeListOfShape MVE;
     GEOMImpl_Block6Explorer::MapShapesAndAncestors
       (aBlockOrComp, TopAbs_VERTEX, TopAbs_EDGE, MVE);
@@ -677,23 +744,11 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetEdge
     return NULL;
   }
 
-  //The GetEdge() doesn't change object so no new function is required.
-  Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
+  Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aResult->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.GetEdge(";
-  TDF_Tool::Entry(theShape->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(thePoint1->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(thePoint2->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  TCollection_AsciiString aNewDescr = aFunction->GetDescription() + "\n";
-  aNewDescr += aDescr;
-  aFunction->SetDescription(aNewDescr);
+  GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetEdge("
+    << theShape << ", " << thePoint1 << ", " << thePoint2 << ")";
 
   SetErrorCode(OK);
   return aResult;
@@ -718,13 +773,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetEdgeNearPoint
 
   TopoDS_Shape aBlockOrComp = theShape->GetValue();
   if (aBlockOrComp.IsNull()) {
-    SetErrorCode("Block or compound is null");
-    return NULL;
-  }
-  if (aBlockOrComp.ShapeType() != TopAbs_SOLID &&
-      aBlockOrComp.ShapeType() != TopAbs_COMPOUND &&
-      aBlockOrComp.ShapeType() != TopAbs_COMPSOLID) {
-    SetErrorCode("Shape is neither a block, nor a compound of blocks");
+    SetErrorCode("Given shape is null");
     return NULL;
   }
 
@@ -740,66 +789,17 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetEdgeNearPoint
 
   //Compute the Edge value
   try {
-    TopoDS_Shape aShape;
-
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
+    TopoDS_Shape aShape = GEOMUtils::GetEdgeNearPoint(aBlockOrComp, aVert);
 
-    // 1. Explode blocks on edges
-    TopTools_MapOfShape mapShape;
-    Standard_Integer nbEdges = 0;
-    TopExp_Explorer exp (aBlockOrComp, TopAbs_EDGE);
-    for (; exp.More(); exp.Next()) {
-      if (mapShape.Add(exp.Current())) {
-        nbEdges++;
-      }
-    }
-
-    mapShape.Clear();
-    Standard_Integer ind = 1;
-    TopTools_Array1OfShape anEdges (1, nbEdges);
-    TColStd_Array1OfReal aDistances (1, nbEdges);
-    for (exp.Init(aBlockOrComp, TopAbs_EDGE); exp.More(); exp.Next()) {
-      if (mapShape.Add(exp.Current())) {
-        TopoDS_Shape anEdge = exp.Current();
-        anEdges(ind) = anEdge;
-
-        // 2. Classify the point relatively each edge
-        BRepExtrema_DistShapeShape aDistTool (aVert, anEdges(ind));
-        if (!aDistTool.IsDone()) {
-          SetErrorCode("Can not find a distance from the given point to one of edges");
-          return NULL;
-        }
-        aDistances(ind) = aDistTool.Value();
-        ind++;
-      }
-    }
-
-    // 3. Define edge, having minimum distance to the point
-    Standard_Real nearest = RealLast(), nbFound = 0;
-    Standard_Real prec = Precision::Confusion();
-    for (ind = 1; ind <= nbEdges; ind++) {
-      if (Abs(aDistances(ind) - nearest) < prec) {
-        nbFound++;
-      } else if (aDistances(ind) < nearest) {
-        nearest = aDistances(ind);
-        aShape = anEdges(ind);
-        nbFound = 1;
-      } else {
-      }
-    }
-    if (nbFound > 1) {
-      SetErrorCode("Multiple edges near the given point are found");
-      return NULL;
-    } else if (nbFound == 0) {
-      SetErrorCode("There are no edges near the given point");
-      return NULL;
-    } else {
-      TopTools_IndexedMapOfShape anIndices;
-      TopExp::MapShapes(aBlockOrComp, anIndices);
-      Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
-      anArray->SetValue(1, anIndices.FindIndex(aShape));
-      aResult = GetEngine()->AddSubShape(theShape, anArray);
-    }
+    TopTools_IndexedMapOfShape anIndices;
+    TopExp::MapShapes(aBlockOrComp, anIndices);
+    Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
+    anArray->SetValue(1, anIndices.FindIndex(aShape));
+    aResult = GetEngine()->AddSubShape(theShape, anArray);
   }
   catch (Standard_Failure) {
     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
@@ -807,21 +807,11 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetEdgeNearPoint
     return NULL;
   }
 
-  //The GetEdgeNearPoint() doesn't change object so no new function is required.
-  Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
+  Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aResult->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.GetEdgeNearPoint(";
-  TDF_Tool::Entry(theShape->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(thePoint->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  TCollection_AsciiString aNewDescr = aFunction->GetDescription() + "\n";
-  aNewDescr += aDescr;
-  aFunction->SetDescription(aNewDescr);
+  GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetEdgeNearPoint("
+                               << theShape << ", " << thePoint << ")";
 
   SetErrorCode(OK);
   return aResult;
@@ -854,12 +844,6 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByPoints
     SetErrorCode("Block or compound is null");
     return NULL;
   }
-  if (aBlockOrComp.ShapeType() != TopAbs_SOLID &&
-      aBlockOrComp.ShapeType() != TopAbs_COMPOUND &&
-      aBlockOrComp.ShapeType() != TopAbs_COMPSOLID) {
-    SetErrorCode("Shape is neither a block, nor a compound of blocks");
-    return NULL;
-  }
 
   TopoDS_Shape anArg1 = thePoint1->GetValue();
   TopoDS_Shape anArg2 = thePoint2->GetValue();
@@ -880,6 +864,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByPoints
 
   //Compute the Face value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     TopoDS_Shape aShape;
 
     TopTools_IndexedDataMapOfShapeListOfShape MVF;
@@ -963,27 +950,12 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByPoints
     return NULL;
   }
 
-  //The GetFaceByPoints() doesn't change object so no new function is required.
-  Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
+  Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aResult->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.GetFaceByPoints(";
-  TDF_Tool::Entry(theShape->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(thePoint1->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(thePoint2->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(thePoint3->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(thePoint4->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  TCollection_AsciiString aNewDescr = aFunction->GetDescription() + "\n";
-  aNewDescr += aDescr;
-  aFunction->SetDescription(aNewDescr);
+  GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceByPoints("
+    << theShape << ", " << thePoint1 << ", " << thePoint2
+      << ", " << thePoint3 << ", " << thePoint4 << ")";
 
   SetErrorCode(OK);
   return aResult;
@@ -1012,12 +984,6 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByEdges
     SetErrorCode("Block or compound is null");
     return NULL;
   }
-  if (aBlockOrComp.ShapeType() != TopAbs_SOLID &&
-      aBlockOrComp.ShapeType() != TopAbs_COMPOUND &&
-      aBlockOrComp.ShapeType() != TopAbs_COMPSOLID) {
-    SetErrorCode("Shape is neither a block, nor a compound of blocks");
-    return NULL;
-  }
 
   TopoDS_Shape anArg1 = theEdge1->GetValue();
   TopoDS_Shape anArg2 = theEdge2->GetValue();
@@ -1033,6 +999,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByEdges
 
   //Compute the Face value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     TopoDS_Shape aShape;
 
     TopTools_IndexedDataMapOfShapeListOfShape MEF;
@@ -1106,23 +1075,11 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByEdges
     return NULL;
   }
 
-  //The GetFaceByEdges() doesn't change object so no new function is required.
-  Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
+  Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aResult->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.GetFaceByEdges(";
-  TDF_Tool::Entry(theShape->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(theEdge1->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(theEdge2->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  TCollection_AsciiString aNewDescr = aFunction->GetDescription() + "\n";
-  aNewDescr += aDescr;
-  aFunction->SetDescription(aNewDescr);
+  GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceByEdges("
+    << theShape << ", " << theEdge1 << ", " << theEdge2 << ")";
 
   SetErrorCode(OK);
   return aResult;
@@ -1167,6 +1124,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetOppositeFace
 
   //Compute the Face value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     TopoDS_Shape aShape;
 
     GEOMImpl_Block6Explorer aBlockTool;
@@ -1185,21 +1145,11 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetOppositeFace
     return NULL;
   }
 
-  //The GetOppositeFace() doesn't change object so no new function is required.
-  Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
+  Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aResult->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.GetOppositeFace(";
-  TDF_Tool::Entry(theShape->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(theFace->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  TCollection_AsciiString aNewDescr = aFunction->GetDescription() + "\n";
-  aNewDescr += aDescr;
-  aFunction->SetDescription(aNewDescr);
+  GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetOppositeFace("
+                               << theShape << ", " << theFace << ")";
 
   SetErrorCode(OK);
   return aResult;
@@ -1227,12 +1177,6 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceNearPoint
     SetErrorCode("Block or compound is null");
     return NULL;
   }
-  if (aBlockOrComp.ShapeType() != TopAbs_SOLID &&
-      aBlockOrComp.ShapeType() != TopAbs_COMPOUND &&
-      aBlockOrComp.ShapeType() != TopAbs_COMPSOLID) {
-    SetErrorCode("Shape is neither a block, nor a compound of blocks");
-    return NULL;
-  }
 
   TopoDS_Shape anArg = thePoint->GetValue();
   if (anArg.IsNull()) {
@@ -1246,55 +1190,67 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceNearPoint
 
   //Compute the Face value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     TopoDS_Shape aShape;
 
     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
     gp_Pnt aPnt = BRep_Tool::Pnt(aVert);
+    Standard_Real PX, PY, PZ;
+    aPnt.Coord(PX, PY, PZ);
 
-    // 1. Explode blocks on faces
-    TopTools_MapOfShape mapShape;
-    Standard_Integer nbFaces = 0;
+    // 1. Classify the point relatively each face
+    Standard_Integer nearest = 2, nbFound = 0;
+    TopTools_DataMapOfShapeInteger mapShapeDist;
     TopExp_Explorer exp (aBlockOrComp, TopAbs_FACE);
     for (; exp.More(); exp.Next()) {
-      if (mapShape.Add(exp.Current())) {
-        nbFaces++;
-      }
-    }
+      TopoDS_Shape aFace = exp.Current();
+
+      if (!mapShapeDist.IsBound(aFace)) {
+        Standard_Integer aDistance = 2;
+
+        // 1.a. Classify relatively Surface
+        Handle(Geom_Surface) aSurf = BRep_Tool::Surface(TopoDS::Face(aFace));
+        Handle(ShapeAnalysis_Surface) aSurfAna = new ShapeAnalysis_Surface (aSurf);
+        gp_Pnt2d p2dOnSurf = aSurfAna->ValueOfUV(aPnt, Precision::Confusion());
+        gp_Pnt p3dOnSurf = aSurfAna->Value(p2dOnSurf);
+        Standard_Real aDist = p3dOnSurf.Distance(aPnt);
+        if (aDist > Precision::Confusion()) {
+          // OUT of Surface
+          aDistance = 1;
+        } else {
+          // 1.b. Classify relatively the face itself
+          BRepClass_FaceClassifier FC (TopoDS::Face(aFace), p2dOnSurf, Precision::Confusion());
+          if (FC.State() == TopAbs_IN) {
+            aDistance = -1;
+          } else if (FC.State() == TopAbs_ON) {
+            aDistance = 0;
+          } else { // OUT
+            aDistance = 1;
+          }
+        }
 
-    mapShape.Clear();
-    Standard_Integer ind = 1;
-    TopTools_Array1OfShape aFaces (1, nbFaces);
-    TColStd_Array1OfInteger aDistances (1, nbFaces);
-    for (exp.Init(aBlockOrComp, TopAbs_FACE); exp.More(); exp.Next()) {
-      if (mapShape.Add(exp.Current())) {
-        TopoDS_Shape aFace = exp.Current();
-        aFaces(ind) = aFace;
-
-        // 2. Classify the point relatively each face
-        BRepClass_FaceClassifier FC (TopoDS::Face(aFace), aPnt, Precision::Confusion());
-        if (FC.State() == TopAbs_IN) {
-          aDistances(ind) = -1;
-        } else if (FC.State() == TopAbs_ON) {
-          aDistances(ind) = 0;
-        } else { // OUT
-          aDistances(ind) = 1;
+        if (aDistance < nearest) {
+          nearest = aDistance;
+          aShape = aFace;
+          nbFound = 1;
+
+          // A first found face, containing the point inside, will be returned.
+          // It is the solution, if there are no
+          // coincident or intersecting faces in the compound.
+          if (nearest == -1) break;
+
+        } else if (aDistance == nearest) {
+          nbFound++;
+        } else {
         }
-        ind++;
-      }
-    }
 
-    // 3. Define face, containing the point or having minimum distance to it
-    Standard_Integer nearest = 2, nbFound = 0;
-    for (ind = 1; ind <= nbFaces; ind++) {
-      if (aDistances(ind) < nearest) {
-        nearest = aDistances(ind);
-        aShape = aFaces(ind);
-        nbFound = 1;
-      } else if (aDistances(ind) == nearest) {
-        nbFound++;
-      } else {
-      }
+        mapShapeDist.Bind(aFace, aDistance);
+      } // if (!mapShapeDist.IsBound(aFace))
     }
+
+    // 2. Define face, containing the point or having minimum distance to it
     if (nbFound > 1) {
       if (nearest == 0) {
         // The point is on boundary of some faces and there are
@@ -1306,40 +1262,63 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceNearPoint
         // The point is outside some faces and there are
         // no faces, having the point inside or on boundary.
         // We will get a nearest face
-        Standard_Real minDist = RealLast();
-        for (ind = 1; ind <= nbFaces; ind++) {
-          if (aDistances(ind) == 1) {
-            BRepExtrema_DistShapeShape aDistTool (aVert, aFaces(ind));
-            if (!aDistTool.IsDone()) {
-              SetErrorCode("Can not find a distance from the given point to one of faces");
-              return NULL;
-            }
-            Standard_Real aDist = aDistTool.Value();
-            if (aDist < minDist) {
-              minDist = aDist;
-              aShape = aFaces(ind);
+        Standard_Real bigReal = RealLast();
+        Standard_Real minDist = bigReal;
+        TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
+        for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
+          if (mapShapeDistIter.Value() == 1) {
+            TopoDS_Shape aFace = mapShapeDistIter.Key();
+            Standard_Real aDist = bigReal;
+
+            // 2.a. Fast check of distance - if point projection on surface is on face
+            Handle(Geom_Surface) aSurf = BRep_Tool::Surface(TopoDS::Face(aFace));
+            Handle(ShapeAnalysis_Surface) aSurfAna = new ShapeAnalysis_Surface (aSurf);
+            gp_Pnt2d p2dOnSurf = aSurfAna->ValueOfUV(aPnt, Precision::Confusion());
+            gp_Pnt p3dOnSurf = aSurfAna->Value(p2dOnSurf);
+            aDist = p3dOnSurf.Distance(aPnt);
+
+            BRepClass_FaceClassifier FC (TopoDS::Face(aFace), p2dOnSurf, Precision::Confusion());
+            if (FC.State() == TopAbs_OUT) {
+              if (aDist < minDist) {
+                // 2.b. Slow check - if point projection on surface is outside of face
+                BRepExtrema_DistShapeShape aDistTool (aVert, aFace);
+                if (!aDistTool.IsDone()) {
+                  SetErrorCode("Can not find a distance from the given point to one of faces");
+                  return NULL;
+                }
+                aDist = aDistTool.Value();
+              } else {
+                aDist = bigReal;
+              }
             }
-          }
-        }
-      } else { // nearest == -1
-        // The point is inside some faces.
-        // We will get a face with nearest center
-        Standard_Real minDist = RealLast();
-        for (ind = 1; ind <= nbFaces; ind++) {
-          if (aDistances(ind) == -1) {
-            GProp_GProps aSystem;
-            BRepGProp::SurfaceProperties(aFaces(ind), aSystem);
-            gp_Pnt aCenterMass = aSystem.CentreOfMass();
 
-            Standard_Real aDist = aCenterMass.Distance(aPnt);
             if (aDist < minDist) {
               minDist = aDist;
-              aShape = aFaces(ind);
+              aShape = aFace;
             }
           }
         }
+      } else { // nearest == -1
+//        // The point is inside some faces.
+//        // We will get a face with nearest center
+//        Standard_Real minDist = RealLast();
+//        TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
+//        for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
+//          if (mapShapeDistIter.Value() == -1) {
+//            TopoDS_Shape aFace = mapShapeDistIter.Key();
+//            GProp_GProps aSystem;
+//            BRepGProp::SurfaceProperties(aFace, aSystem);
+//            gp_Pnt aCenterMass = aSystem.CentreOfMass();
+//
+//            Standard_Real aDist = aCenterMass.Distance(aPnt);
+//            if (aDist < minDist) {
+//              minDist = aDist;
+//              aShape = aFace;
+//            }
+//          }
+//        }
       }
-    }
+    } // if (nbFound > 1)
 
     if (nbFound == 0) {
       SetErrorCode("There are no faces near the given point");
@@ -1358,21 +1337,11 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceNearPoint
     return NULL;
   }
 
-  //The GetFaceNearPoint() doesn't change object so no new function is required.
-  Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
+  Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aResult->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.GetFaceNearPoint(";
-  TDF_Tool::Entry(theShape->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(thePoint->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  TCollection_AsciiString aNewDescr = aFunction->GetDescription() + "\n";
-  aNewDescr += aDescr;
-  aFunction->SetDescription(aNewDescr);
+  GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceNearPoint("
+                               << theShape << ", " << thePoint << ")";
 
   SetErrorCode(OK);
   return aResult;
@@ -1400,12 +1369,6 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByNormale
     SetErrorCode("Block or compound is null");
     return NULL;
   }
-  if (aBlockOrComp.ShapeType() != TopAbs_SOLID &&
-      aBlockOrComp.ShapeType() != TopAbs_COMPOUND &&
-      aBlockOrComp.ShapeType() != TopAbs_COMPSOLID) {
-    SetErrorCode("Shape is neither a block, nor a compound of blocks");
-    return NULL;
-  }
 
   TopoDS_Shape anArg = theVector->GetValue();
   if (anArg.IsNull()) {
@@ -1419,6 +1382,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByNormale
 
   //Compute the Face value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     TopoDS_Shape aShape;
 
     TopoDS_Edge anEdge = TopoDS::Edge(anArg);
@@ -1501,21 +1467,11 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByNormale
     return NULL;
   }
 
-  //The GetFaceByNormale() doesn't change object so no new function is required.
-  Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
+  Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aResult->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.GetFaceByNormale(";
-  TDF_Tool::Entry(theShape->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(theVector->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  TCollection_AsciiString aNewDescr = aFunction->GetDescription() + "\n";
-  aNewDescr += aDescr;
-  aFunction->SetDescription(aNewDescr);
+  GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceByNormale("
+    << theShape << ", " << theVector << ")";
 
   SetErrorCode(OK);
   return aResult;
@@ -1523,54 +1479,188 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByNormale
 
 //=============================================================================
 /*!
- *  IsCompoundOfBlocks
+ *  GetShapesNearPoint
  */
 //=============================================================================
-Standard_Boolean GEOMImpl_IBlocksOperations::IsCompoundOfBlocks
-                                                (Handle(GEOM_Object)    theCompound,
-                                                 const Standard_Integer theMinNbFaces,
-                                                 const Standard_Integer theMaxNbFaces,
-                                                 Standard_Integer&      theNbBlocks)
+Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetShapesNearPoint
+                                         (Handle(GEOM_Object)    theShape,
+                                          Handle(GEOM_Object)    thePoint,
+                                          const Standard_Integer theShapeType,
+                                          const Standard_Real    theConstTolerance)
 {
   SetErrorCode(KO);
-  Standard_Boolean isCompOfBlocks = Standard_False;
-  theNbBlocks = 0;
 
-  if (theCompound.IsNull()) return isCompOfBlocks;
-  TopoDS_Shape aBlockOrComp = theCompound->GetValue();
+  // New object
+  Handle(GEOM_Object) aResult;
 
-  //Check
-  isCompOfBlocks = Standard_True;
+  // Arguments
+  if (theShape.IsNull() || thePoint.IsNull()) return NULL;
+
+  TopoDS_Shape aBlockOrComp = theShape->GetValue();
+  if (aBlockOrComp.IsNull()) {
+    SetErrorCode("Block or compound is null");
+    return NULL;
+  }
+
+  TopoDS_Shape anArg = thePoint->GetValue();
+  if (anArg.IsNull()) {
+    SetErrorCode("Null shape is given as argument");
+    return NULL;
+  }
+  if (anArg.ShapeType() != TopAbs_VERTEX) {
+    SetErrorCode("Element for face identification is not a vertex");
+    return NULL;
+  }
+
+  if (theShapeType < TopAbs_SOLID || TopAbs_VERTEX < theShapeType) {
+    SetErrorCode("Invalid type of result is requested");
+    return NULL;
+  }
+  
+  Standard_Real theTolerance = theConstTolerance;
+  if (theTolerance < Precision::Confusion()) {
+    theTolerance = Precision::Confusion();
+  }
+
+  // Compute the result
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
+    TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
+
     TopTools_MapOfShape mapShape;
-    TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
+    Standard_Integer nbEdges = 0;
+    TopExp_Explorer exp (aBlockOrComp, TopAbs_ShapeEnum(theShapeType));
     for (; exp.More(); exp.Next()) {
       if (mapShape.Add(exp.Current())) {
-        TopoDS_Shape aSolid = exp.Current();
+        nbEdges++;
+      }
+    }
 
-        TopTools_MapOfShape mapFaces;
-        TopExp_Explorer expF (aSolid, TopAbs_FACE);
-        Standard_Integer nbFaces = 0;
-        for (; expF.More(); expF.Next()) {
-          if (mapFaces.Add(expF.Current())) {
-            nbFaces++;
-            if (nbFaces > theMaxNbFaces) {
-              isCompOfBlocks = Standard_False;
-              break;
-            }
-          }
+    if (nbEdges == 0) {
+      SetErrorCode("Given shape contains no sub-shapes of requested type");
+      return NULL;
+    }
+
+    // Calculate distances and find min
+    mapShape.Clear();
+    Standard_Integer ind = 1;
+    Standard_Real aMinDist = RealLast();
+    TopTools_Array1OfShape anEdges (1, nbEdges);
+    TColStd_Array1OfReal aDistances (1, nbEdges);
+    for (exp.Init(aBlockOrComp, TopAbs_ShapeEnum(theShapeType)); exp.More(); exp.Next()) {
+      if (mapShape.Add(exp.Current())) {
+        TopoDS_Shape anEdge = exp.Current();
+        anEdges(ind) = anEdge;
+
+        BRepExtrema_DistShapeShape aDistTool (aVert, anEdges(ind));
+        if (!aDistTool.IsDone()) {
+          SetErrorCode("Can not find a distance from the given point to one of sub-shapes");
+          return NULL;
         }
-        if (nbFaces < theMinNbFaces || theMaxNbFaces < nbFaces) {
-          isCompOfBlocks = Standard_False;
-        } else {
-          theNbBlocks++;
+        aDistances(ind) = aDistTool.Value();
+        if (aDistances(ind) < aMinDist) {
+          aMinDist = aDistances(ind);
         }
+        ind++;
       }
     }
-  }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+
+    if (aMinDist < RealLast()) {
+      // Collect sub-shapes with distance < (aMinDist + theTolerance)
+      int nbSubShapes = 0;
+      TopTools_Array1OfShape aNearShapes (1, nbEdges);
+      for (ind = 1; ind <= nbEdges; ind++) {
+        if (aDistances(ind) < aMinDist + theTolerance) {
+          nbSubShapes++;
+          aNearShapes(nbSubShapes) = anEdges(ind);
+        }
+      }
+
+      // Add sub-shape
+      TopTools_IndexedMapOfShape anIndices;
+      TopExp::MapShapes(aBlockOrComp, anIndices);
+      Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger (1, nbSubShapes);
+      for (ind = 1; ind <= nbSubShapes; ind++) {
+        anArray->SetValue(ind, anIndices.FindIndex(aNearShapes(ind)));
+      }
+      aResult = GetEngine()->AddSubShape(theShape, anArray);
+    }
+  }
+  catch (Standard_Failure) {
+    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
+    SetErrorCode(aFail->GetMessageString());
+    return NULL;
+  }
+
+  if (aResult.IsNull())
+    return NULL;
+
+  Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
+
+  //Make a Python command
+  GEOM::TPythonDump(aFunction)
+    << aResult << " = geompy.GetShapesNearPoint(" << theShape << ", " << thePoint
+    << ", " << TopAbs_ShapeEnum(theShapeType) << ", " << theTolerance << ")";
+
+  SetErrorCode(OK);
+  return aResult;
+}
+
+//=============================================================================
+/*!
+ *  IsCompoundOfBlocks
+ */
+//=============================================================================
+Standard_Boolean GEOMImpl_IBlocksOperations::IsCompoundOfBlocks
+                                                (Handle(GEOM_Object)    theCompound,
+                                                 const Standard_Integer theMinNbFaces,
+                                                 const Standard_Integer theMaxNbFaces,
+                                                 Standard_Integer&      theNbBlocks)
+{
+  SetErrorCode(KO);
+  Standard_Boolean isCompOfBlocks = Standard_False;
+  theNbBlocks = 0;
+
+  if (theCompound.IsNull()) return isCompOfBlocks;
+  TopoDS_Shape aBlockOrComp = theCompound->GetValue();
+
+  //Check
+  isCompOfBlocks = Standard_True;
+  try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
+    TopTools_MapOfShape mapShape;
+    TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
+    for (; exp.More(); exp.Next()) {
+      if (mapShape.Add(exp.Current())) {
+        TopoDS_Shape aSolid = exp.Current();
+
+        TopTools_MapOfShape mapFaces;
+        TopExp_Explorer expF (aSolid, TopAbs_FACE);
+        Standard_Integer nbFaces = 0;
+        for (; expF.More(); expF.Next()) {
+          if (mapFaces.Add(expF.Current())) {
+            nbFaces++;
+            if (nbFaces > theMaxNbFaces) {
+              isCompOfBlocks = Standard_False;
+              break;
+            }
+          }
+        }
+        if (nbFaces < theMinNbFaces || theMaxNbFaces < nbFaces) {
+          isCompOfBlocks = Standard_False;
+        } else {
+          theNbBlocks++;
+        }
+      }
+    }
+  }
+  catch (Standard_Failure) {
+    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
+    SetErrorCode(aFail->GetMessageString());
     return isCompOfBlocks;
   }
 
@@ -1583,9 +1673,139 @@ Standard_Boolean GEOMImpl_IBlocksOperations::IsCompoundOfBlocks
  *  Set of functions, used by CheckCompoundOfBlocks() method
  */
 //=============================================================================
-void AddBlocksFrom (const TopoDS_Shape&  theShape,
-                    TopTools_ListOfShape& BLO,
-                    TopTools_ListOfShape& NOT)
+void GEOMImpl_IBlocksOperations::AddBlocksFrom (const TopoDS_Shape&   theShape,
+                                                TopTools_ListOfShape& BLO,
+                                                TopTools_ListOfShape& NOT,
+                                                TopTools_ListOfShape& EXT,
+                                                TopTools_ListOfShape& NOQ)
+{
+  TopAbs_ShapeEnum aType = theShape.ShapeType();
+  switch (aType) {
+  case TopAbs_COMPOUND:
+  case TopAbs_COMPSOLID:
+    {
+      TopoDS_Iterator It (theShape);
+      for (; It.More(); It.Next()) {
+        AddBlocksFrom(It.Value(), BLO, NOT, EXT, NOQ);
+      }
+    }
+    break;
+  case TopAbs_SOLID:
+    {
+      // Check, if there are seam or degenerated edges
+      BlockFix_CheckTool aTool;
+      aTool.SetShape(theShape);
+      aTool.Perform();
+      if (aTool.NbPossibleBlocks() > 0) {
+        EXT.Append(theShape);
+      } else {
+        // Count faces and edges in each face to recognize blocks
+        TopTools_MapOfShape mapFaces;
+        Standard_Integer nbFaces = 0;
+        Standard_Boolean hasNonQuadr = Standard_False;
+        TopExp_Explorer expF (theShape, TopAbs_FACE);
+
+        for (; expF.More(); expF.Next()) {
+          if (mapFaces.Add(expF.Current())) {
+            nbFaces++;
+            //0021483//if (nbFaces > 6) break;
+
+            // get wire
+            TopoDS_Shape aF = expF.Current();
+            TopExp_Explorer wires (aF, TopAbs_WIRE);
+            if (!wires.More()) {
+              // no wire in the face
+              hasNonQuadr = Standard_True;
+              NOQ.Append(aF);//0021483
+              //0021483//break;
+              continue;
+            }
+            TopoDS_Shape aWire = wires.Current();
+            wires.Next();
+            if (wires.More()) {
+              // multiple wires in the face
+              hasNonQuadr = Standard_True;
+              NOQ.Append(aF);//0021483
+              //0021483//break;
+              continue;
+            }
+
+            // Check number of edges in the face
+            Standard_Integer nbEdges = 0;
+            TopTools_MapOfShape mapEdges;
+            TopExp_Explorer expW (aWire, TopAbs_EDGE);
+            for (; expW.More(); expW.Next()) {
+              if (mapEdges.Add(expW.Current())) {
+                nbEdges++;
+                if (nbEdges > 4) break;
+              }
+            }
+            if (nbEdges != 4) {
+              hasNonQuadr = Standard_True;
+              NOQ.Append(aF);//0021483
+            }
+          }
+        }
+
+        if (nbFaces == 6 && !hasNonQuadr) {
+          BLO.Append(theShape);
+        } else {
+          NOT.Append(theShape);
+        }
+      }
+    }
+    break;
+  case TopAbs_SHELL: //0021483
+  case TopAbs_FACE: //0021483
+    {
+      // Count edges in each face
+      TopTools_MapOfShape mapFaces;
+      TopExp_Explorer expF (theShape, TopAbs_FACE);
+      for (; expF.More(); expF.Next()) {
+        if (mapFaces.Add(expF.Current())) {
+          // get wire
+          TopoDS_Shape aF = expF.Current();
+          TopExp_Explorer wires (aF, TopAbs_WIRE);
+          if (!wires.More()) {
+            // no wire in the face
+            NOQ.Append(aF);//0021483
+            continue;
+          }
+          TopoDS_Shape aWire = wires.Current();
+          wires.Next();
+          if (wires.More()) {
+            // multiple wires in the face
+            NOQ.Append(aF);//0021483
+            continue;
+          }
+
+          // Check number of edges in the face
+          Standard_Integer nbEdges = 0;
+          TopTools_MapOfShape mapEdges;
+          TopExp_Explorer expW (aWire, TopAbs_EDGE);
+          for (; expW.More(); expW.Next()) {
+            if (mapEdges.Add(expW.Current())) {
+              nbEdges++;
+              if (nbEdges > 4) break;
+            }
+          }
+          if (nbEdges != 4) {
+            NOQ.Append(aF);//0021483
+          }
+        }
+      }
+    }
+    break;
+  default:
+    NOT.Append(theShape);
+  }
+}
+
+void AddBlocksFromOld (const TopoDS_Shape&   theShape,
+                       TopTools_ListOfShape& BLO,
+                       TopTools_ListOfShape& NOT,
+                       TopTools_ListOfShape& DEG,
+                       TopTools_ListOfShape& SEA)
 {
   TopAbs_ShapeEnum aType = theShape.ShapeType();
   switch (aType) {
@@ -1594,7 +1814,7 @@ void AddBlocksFrom (const TopoDS_Shape&  theShape,
     {
       TopoDS_Iterator It (theShape);
       for (; It.More(); It.Next()) {
-        AddBlocksFrom(It.Value(), BLO, NOT);
+        AddBlocksFromOld(It.Value(), BLO, NOT, DEG, SEA);
       }
     }
     break;
@@ -1603,25 +1823,67 @@ void AddBlocksFrom (const TopoDS_Shape&  theShape,
       TopTools_MapOfShape mapFaces;
       TopExp_Explorer expF (theShape, TopAbs_FACE);
       Standard_Integer nbFaces = 0;
-      Standard_Integer nbEdges = 0;
+      Standard_Boolean hasNonQuadr = Standard_False;
+      Standard_Boolean hasDegenerated = Standard_False;
+      Standard_Boolean hasSeam = Standard_False;
       for (; expF.More(); expF.Next()) {
         if (mapFaces.Add(expF.Current())) {
           nbFaces++;
           if (nbFaces > 6) break;
 
           // Check number of edges in the face
+          Standard_Integer nbEdges = 0;
+          TopTools_MapOfShape mapEdges;
+
+          // get wire
           TopoDS_Shape aF = expF.Current();
-          TopExp_Explorer expE (aF, TopAbs_EDGE);
-          nbEdges = 0;
-          for (; expE.More(); expE.Next()) {
-            nbEdges++;
-            if (nbEdges > 4) break;
+          TopExp_Explorer wires (aF, TopAbs_WIRE);
+          if (!wires.More()) {
+            // no wire in the face
+            hasNonQuadr = Standard_True;
+            break;
+          }
+          TopoDS_Shape aWire = wires.Current();
+          wires.Next();
+          if (wires.More()) {
+            // multiple wires in the face
+            hasNonQuadr = Standard_True;
+            break;
+          }
+
+          // iterate on wire
+          BRepTools_WireExplorer aWE (TopoDS::Wire(aWire), TopoDS::Face(aF));
+          for (; aWE.More(); aWE.Next(), nbEdges++) {
+            if (BRep_Tool::Degenerated(aWE.Current())) {
+              // degenerated edge found
+              hasDegenerated = Standard_True;
+//              break;
+            }
+            if (mapEdges.Contains(aWE.Current())) {
+              // seam edge found
+              hasSeam = Standard_True;
+//              break;
+            }
+            mapEdges.Add(aWE.Current());
+          }
+          if (nbEdges != 4) {
+            hasNonQuadr = Standard_True;
           }
-          if (nbEdges != 4) break;
         }
       }
-      if (nbFaces == 6 && nbEdges == 4) {
-        BLO.Append(theShape);
+      if (nbFaces == 6) {
+        if (hasDegenerated || hasSeam) {
+          if (hasDegenerated) {
+            DEG.Append(theShape);
+          }
+          if (hasSeam) {
+            SEA.Append(theShape);
+          }
+        } else if (hasNonQuadr) {
+          NOT.Append(theShape);
+        } else {
+          BLO.Append(theShape);
+        }
       } else {
         NOT.Append(theShape);
       }
@@ -1649,20 +1911,13 @@ Standard_Integer BlocksRelation (const TopoDS_Shape& theBlock1,
   Bnd_Box B1, B2;
   BRepBndLib::Add(theBlock1, B1);
   BRepBndLib::Add(theBlock2, B2);
-//  BRepBndLib::AddClose(theBlock1, B1);
-//  BRepBndLib::AddClose(theBlock2, B2);
   B1.Get(Xmin1, Ymin1, Zmin1, Xmax1, Ymax1, Zmax1);
   B2.Get(Xmin2, Ymin2, Zmin2, Xmax2, Ymax2, Zmax2);
   if (Xmax2 < Xmin1 || Xmax1 < Xmin2 ||
       Ymax2 < Ymin1 || Ymax1 < Ymin2 ||
       Zmax2 < Zmin1 || Zmax1 < Zmin2) {
-//  Standard_Real prec = Precision::Confusion();
-//  if (prec < Xmin1 - Xmax2 || prec < Xmin2 - Xmax1 ||
-//      prec < Ymin1 - Ymax2 || prec < Ymin2 - Ymax1 ||
-//      prec < Zmin1 - Zmax2 || prec < Zmin2 - Zmax1) {
     return REL_NOT_CONNECTED;
   }
-  // to be done
 
   BRepExtrema_DistShapeShape dst (theBlock1, theBlock2);
   if (!dst.IsDone()) {
@@ -1863,12 +2118,12 @@ Standard_Boolean HasAnyConnection (const Standard_Integer         theBlockIndex,
 
 //=============================================================================
 /*!
- *  CheckCompoundOfBlocks
+ *  CheckCompoundOfBlocksOld
  */
 //=============================================================================
-Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
+Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocksOld
                                                 (Handle(GEOM_Object) theCompound,
-                                                 list<BCError>&      theErrors)
+                                                 std::list<BCError>&      theErrors)
 {
   SetErrorCode(KO);
 
@@ -1883,20 +2138,40 @@ Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
 
   // 1. Report non-blocks
   TopTools_ListOfShape NOT; // Not blocks
+  TopTools_ListOfShape DEG; // Hexahedral solids, having degenerated edges
+  TopTools_ListOfShape SEA; // Hexahedral solids, having seam edges
   TopTools_ListOfShape BLO; // All blocks from the given compound
-  AddBlocksFrom(aBlockOrComp, BLO, NOT);
+  AddBlocksFromOld(aBlockOrComp, BLO, NOT, DEG, SEA);
 
   if (NOT.Extent() > 0) {
     isCompOfBlocks = Standard_False;
     BCError anErr;
     anErr.error = NOT_BLOCK;
-    TopTools_ListIteratorOfListOfShape NOTit (NOT);
-    for (; NOTit.More(); NOTit.Next()) {
-      anErr.incriminated.push_back(anIndices.FindIndex(NOTit.Value()));
+    TopTools_ListIteratorOfListOfShape it (NOT);
+    for (; it.More(); it.Next()) {
+      anErr.incriminated.push_back(anIndices.FindIndex(it.Value()));
     }
     theErrors.push_back(anErr);
   }
 
+  if (DEG.Extent() > 0 || SEA.Extent() > 0) {
+    isCompOfBlocks = Standard_False;
+    BCError anErr;
+    anErr.error = EXTRA_EDGE;
+
+    TopTools_ListIteratorOfListOfShape itDEG (DEG);
+    for (; itDEG.More(); itDEG.Next()) {
+      anErr.incriminated.push_back(anIndices.FindIndex(itDEG.Value()));
+    }
+
+    TopTools_ListIteratorOfListOfShape itSEA (SEA);
+    for (; itSEA.More(); itSEA.Next()) {
+      anErr.incriminated.push_back(anIndices.FindIndex(itSEA.Value()));
+    }
+
+    theErrors.push_back(anErr);
+  }
+
   Standard_Integer nbBlocks = BLO.Extent();
   if (nbBlocks == 0) {
     isCompOfBlocks = Standard_False;
@@ -1962,6 +2237,7 @@ Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
   TColStd_MapOfInteger aCurrentSet;
   for (ibl = 1; ibl <= nbBlocks; ibl++) {
     if (!aProcessedMap.Contains(ibl)) {
+      aCurrentSet.Clear();
       FindConnected(ibl, aRelations, aProcessedMap, aCurrentSet);
       if (aCurrentSet.Extent() > aLargestSet.Extent()) {
         aLargestSet = aCurrentSet;
@@ -1999,34 +2275,37 @@ Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
 //=============================================================================
 TCollection_AsciiString GEOMImpl_IBlocksOperations::PrintBCErrors
                                                 (Handle(GEOM_Object)  theCompound,
-                                                 const list<BCError>& theErrors)
+                                                 const std::list<BCError>& theErrors)
 {
   TCollection_AsciiString aDescr;
 
-  list<BCError>::const_iterator errIt = theErrors.begin();
+  std::list<BCError>::const_iterator errIt = theErrors.begin();
   int i = 0;
   for (; errIt != theErrors.end(); i++, errIt++) {
     BCError errStruct = *errIt;
 
     switch (errStruct.error) {
     case NOT_BLOCK:
-      aDescr += "\nNot a Blocks: ";
+      aDescr += "\n\tNot a Blocks: ";
+      break;
+    case EXTRA_EDGE:
+      aDescr += "\n\tHexahedral solids with degenerated and/or seam edges: ";
       break;
     case INVALID_CONNECTION:
-      aDescr += "\nInvalid connection between two blocks: ";
+      aDescr += "\n\tInvalid connection between two blocks: ";
       break;
     case NOT_CONNECTED:
-      aDescr += "\nBlocks, not connected with main body: ";
+      aDescr += "\n\tBlocks, not connected with main body: ";
       break;
     case NOT_GLUED:
-      aDescr += "\nNot glued blocks: ";
+      aDescr += "\n\tNot glued blocks: ";
       break;
     default:
       break;
     }
 
-    list<int> sshList = errStruct.incriminated;
-    list<int>::iterator sshIt = sshList.begin();
+    std::list<int> sshList = errStruct.incriminated;
+    std::list<int>::iterator sshIt = sshList.begin();
     int jj = 0;
     for (; sshIt != sshList.end(); jj++, sshIt++) {
       if (jj > 0)
@@ -2038,6 +2317,448 @@ TCollection_AsciiString GEOMImpl_IBlocksOperations::PrintBCErrors
   return aDescr;
 }
 
+//=============================================================================
+/*!
+ *  CheckCompoundOfBlocks
+ */
+//=============================================================================
+Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
+                                              (Handle(GEOM_Object) theCompound,
+                                               std::list<BCError>& theErrors)
+{
+  SetErrorCode(KO);
+
+  if (theCompound.IsNull()) return Standard_False;
+  TopoDS_Shape aBlockOrComp = theCompound->GetValue();
+
+  Standard_Boolean isCompOfBlocks = Standard_True;
+
+  // Map sub-shapes and their indices
+  TopTools_IndexedMapOfShape anIndices;
+  TopExp::MapShapes(aBlockOrComp, anIndices);
+
+  // 1. Separate blocks from non-blocks
+  TopTools_ListOfShape NOT; // Not blocks
+  TopTools_ListOfShape EXT; // Hexahedral solids, having degenerated and/or seam edges
+  TopTools_ListOfShape BLO; // All blocks from the given compound
+  TopTools_ListOfShape NOQ; // All non-quadrangular faces
+  AddBlocksFrom(aBlockOrComp, BLO, NOT, EXT, NOQ);
+
+  // Report non-blocks
+  if (NOT.Extent() > 0) {
+    isCompOfBlocks = Standard_False;
+    BCError anErr;
+    anErr.error = NOT_BLOCK;
+    TopTools_ListIteratorOfListOfShape it (NOT);
+    for (; it.More(); it.Next()) {
+      anErr.incriminated.push_back(anIndices.FindIndex(it.Value()));
+    }
+    theErrors.push_back(anErr);
+  }
+
+  // Report solids, having degenerated and/or seam edges
+  if (EXT.Extent() > 0) {
+    isCompOfBlocks = Standard_False;
+    BCError anErr;
+    anErr.error = EXTRA_EDGE;
+    TopTools_ListIteratorOfListOfShape it (EXT);
+    for (; it.More(); it.Next()) {
+      anErr.incriminated.push_back(anIndices.FindIndex(it.Value()));
+    }
+    theErrors.push_back(anErr);
+  }
+
+  Standard_Integer nbBlocks = BLO.Extent();
+  if (nbBlocks == 0) {
+    isCompOfBlocks = Standard_False;
+    SetErrorCode(OK);
+    return isCompOfBlocks;
+  }
+  if (nbBlocks == 1) {
+    SetErrorCode(OK);
+    return isCompOfBlocks;
+  }
+
+  // Prepare data for 2. and 3.
+  TColStd_Array2OfInteger aRelations (1, nbBlocks, 1, nbBlocks);
+  aRelations.Init(REL_NOT_CONNECTED);
+
+  TopTools_IndexedMapOfShape mapBlocks;
+
+  BRep_Builder BB;
+  TopoDS_Compound aComp;
+  BB.MakeCompound(aComp);
+
+  TopTools_ListIteratorOfListOfShape BLOit (BLO);
+  for (; BLOit.More(); BLOit.Next()) {
+    mapBlocks.Add(BLOit.Value());
+    BB.Add(aComp, BLOit.Value());
+  }
+
+  // 2. Find glued blocks (having shared faces)
+  TopTools_IndexedDataMapOfShapeListOfShape mapFaceBlocks;
+  GEOMImpl_Block6Explorer::MapShapesAndAncestors
+    (aComp, TopAbs_FACE, TopAbs_SOLID, mapFaceBlocks);
+
+  Standard_Integer prevInd = 0, curInd = 0;
+  Standard_Integer ind = 1, nbFaces = mapFaceBlocks.Extent();
+  for (; ind <= nbFaces; ind++) {
+    const TopTools_ListOfShape& aGluedBlocks = mapFaceBlocks.FindFromIndex(ind);
+    if (aGluedBlocks.Extent() > 1) { // Shared face found
+      TopTools_ListIteratorOfListOfShape aGluedBlocksIt (aGluedBlocks);
+      TopoDS_Shape prevBlock, curBlock;
+      for (; aGluedBlocksIt.More(); aGluedBlocksIt.Next()) {
+        curBlock = aGluedBlocksIt.Value();
+        if (!prevBlock.IsNull()) {
+          prevInd = mapBlocks.FindIndex(prevBlock);
+          curInd  = mapBlocks.FindIndex(curBlock);
+          aRelations.SetValue(prevInd, curInd, REL_OK);
+          aRelations.SetValue(curInd, prevInd, REL_OK);
+        }
+        prevBlock = curBlock;
+      }
+    }
+  }
+
+  // 3. Find not glued blocks
+  GEOMAlgo_GlueAnalyser aGD;
+
+  aGD.SetShape(aComp);
+  aGD.SetTolerance(Precision::Confusion());
+  aGD.SetCheckGeometry(Standard_True);
+  aGD.Perform();
+
+  Standard_Integer iErr, iWrn;
+  iErr = aGD.ErrorStatus();
+  if (iErr) {
+    SetErrorCode("Error in GEOMAlgo_GlueAnalyser");
+    return isCompOfBlocks;
+  }
+  iWrn = aGD.WarningStatus();
+  if (iWrn) {
+    MESSAGE("Warning in GEOMAlgo_GlueAnalyser");
+  }
+
+  // Report not glued blocks
+  if (aGD.HasSolidsToGlue()) {
+    isCompOfBlocks = Standard_False;
+    Standard_Integer aSx1Ind, aSx2Ind;
+
+    const GEOMAlgo_ListOfCoupleOfShapes& aLCS = aGD.SolidsToGlue();
+    GEOMAlgo_ListIteratorOfListOfCoupleOfShapes aItCS (aLCS);
+    for (; aItCS.More(); aItCS.Next()) {
+      const GEOMAlgo_CoupleOfShapes& aCS = aItCS.Value();
+      const TopoDS_Shape& aSx1 = aCS.Shape1();
+      const TopoDS_Shape& aSx2 = aCS.Shape2();
+
+      aSx1Ind = mapBlocks.FindIndex(aSx1);
+      aSx2Ind = mapBlocks.FindIndex(aSx2);
+      aRelations.SetValue(aSx1Ind, aSx2Ind, NOT_GLUED);
+      aRelations.SetValue(aSx2Ind, aSx1Ind, NOT_GLUED);
+
+      BCError anErr;
+      anErr.error = NOT_GLUED;
+      anErr.incriminated.push_back(anIndices.FindIndex(aSx1));
+      anErr.incriminated.push_back(anIndices.FindIndex(aSx2));
+      theErrors.push_back(anErr);
+    }
+  }
+
+  // 4. Find largest set of connected (good connection or not glued) blocks
+  Standard_Integer ibl = 1;
+  TColStd_MapOfInteger aProcessedMap;
+  TColStd_MapOfInteger aLargestSet;
+  TColStd_MapOfInteger aCurrentSet;
+  for (ibl = 1; ibl <= nbBlocks; ibl++) {
+    if (!aProcessedMap.Contains(ibl)) {
+      aCurrentSet.Clear();
+      FindConnected(ibl, aRelations, aProcessedMap, aCurrentSet);
+      if (aCurrentSet.Extent() > aLargestSet.Extent()) {
+        aLargestSet = aCurrentSet;
+      }
+    }
+  }
+
+  // 5. Report all blocks, isolated from <aLargestSet>
+  BCError anErr;
+  anErr.error = NOT_CONNECTED;
+  Standard_Boolean hasIsolated = Standard_False;
+  for (ibl = 1; ibl <= nbBlocks; ibl++) {
+    if (!aLargestSet.Contains(ibl)) {
+      aProcessedMap.Clear();
+      if (!HasAnyConnection(ibl, aLargestSet, aRelations, aProcessedMap)) {
+        // report connection absence
+        hasIsolated = Standard_True;
+        anErr.incriminated.push_back(anIndices.FindIndex(mapBlocks.FindKey(ibl)));
+      }
+    }
+  }
+  if (hasIsolated) {
+    isCompOfBlocks = Standard_False;
+    theErrors.push_back(anErr);
+  }
+
+  SetErrorCode(OK);
+  return isCompOfBlocks;
+}
+
+//=============================================================================
+/*!
+ *  GetNonBlocks
+ */
+//=============================================================================
+Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetNonBlocks
+                                     (Handle(GEOM_Object) theShape,
+                                      Handle(GEOM_Object)& theNonQuads)
+{
+  SetErrorCode(KO);
+
+  if (theShape.IsNull()) return NULL;
+  TopoDS_Shape aShape = theShape->GetValue();
+
+  // Separate blocks from non-blocks
+  TopTools_ListOfShape BLO; // All blocks from the given compound
+  TopTools_ListOfShape NOT; // Not blocks
+  TopTools_ListOfShape EXT; // Hexahedral solids, having degenerated and/or seam edges
+  TopTools_ListOfShape NOQ; // All non-quadrangular faces
+  AddBlocksFrom(aShape, BLO, NOT, EXT, NOQ);
+
+  if (NOT.IsEmpty() && EXT.IsEmpty() && NOQ.IsEmpty()) {
+    SetErrorCode("NOT_FOUND_ANY");
+    return NULL;
+  }
+
+  // Map sub-shapes and their indices
+  TopTools_IndexedMapOfShape anIndices;
+  TopExp::MapShapes(aShape, anIndices);
+
+  // Non-blocks
+  Handle(GEOM_Object) aNonBlocks;
+  if (NOT.Extent() > 0 || EXT.Extent() > 0) {
+    Handle(TColStd_HArray1OfInteger) anArray =
+      new TColStd_HArray1OfInteger (1, NOT.Extent() + EXT.Extent());
+    Standard_Integer ii = 1;
+    TopTools_ListIteratorOfListOfShape it1 (NOT);
+    for (; it1.More(); it1.Next(), ii++) {
+      anArray->SetValue(ii, anIndices.FindIndex(it1.Value()));
+    }
+    TopTools_ListIteratorOfListOfShape it2 (EXT);
+    for (; it2.More(); it2.Next(), ii++) {
+      anArray->SetValue(ii, anIndices.FindIndex(it2.Value()));
+    }
+
+    aNonBlocks = GetEngine()->AddSubShape(theShape, anArray);
+    if (aNonBlocks.IsNull()) {
+      SetErrorCode("Error in algorithm: result found, but cannot be returned.");
+      return NULL;
+    }
+    aNonBlocks->SetType(GEOM_GROUP);
+    TDF_Label aFreeLabel = aNonBlocks->GetFreeLabel();
+    TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)TopAbs_SOLID);
+  }
+
+  // Non-quadrangles
+  if (NOQ.Extent() > 0) {
+    Handle(TColStd_HArray1OfInteger) anArray =
+      new TColStd_HArray1OfInteger (1, NOQ.Extent());
+    Standard_Integer ii = 1;
+    TopTools_ListIteratorOfListOfShape it1 (NOQ);
+    for (; it1.More(); it1.Next(), ii++) {
+      anArray->SetValue(ii, anIndices.FindIndex(it1.Value()));
+    }
+
+    theNonQuads = GetEngine()->AddSubShape(theShape, anArray);
+    if (theNonQuads.IsNull()) {
+      SetErrorCode("Error in algorithm: result found, but cannot be returned.");
+      return NULL;
+    }
+    theNonQuads->SetType(GEOM_GROUP);
+    TDF_Label aFreeLabel = theNonQuads->GetFreeLabel();
+    TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)TopAbs_FACE);
+  }
+
+  //Make a Python command
+  Handle(GEOM_Function) aMainShape = theShape->GetLastFunction();
+  GEOM::TPythonDump pd (aMainShape, /*append=*/true);
+  pd << "(";
+  if (aNonBlocks.IsNull())
+    pd << "no_bad_solids";
+  else
+    pd << aNonBlocks;
+  pd << ", ";
+  if (theNonQuads.IsNull())
+    pd << "no_bad_faces";
+  else
+    pd << theNonQuads;
+  pd << ") = geompy.GetNonBlocks(" << theShape << ")";
+
+  SetErrorCode(OK);
+  return aNonBlocks;
+}
+
+//=============================================================================
+/*!
+ *  RemoveExtraEdges
+ */
+//=============================================================================
+Handle(GEOM_Object) GEOMImpl_IBlocksOperations::RemoveExtraEdges
+                                     (Handle(GEOM_Object) theObject,
+                                      const Standard_Integer theOptimumNbFaces)
+{
+  SetErrorCode(KO);
+
+  if (theObject.IsNull()) return NULL;
+
+  Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
+  if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be fixed
+
+  //Add a new Copy object
+  Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
+
+  //Add a function
+  Handle(GEOM_Function) aFunction =
+    aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_REMOVE_EXTRA);
+
+  //Check if the function is set correctly
+  if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
+
+  GEOMImpl_IBlockTrsf aTI (aFunction);
+  aTI.SetOriginal(aLastFunction);
+  aTI.SetOptimumNbFaces(theOptimumNbFaces);
+
+  //Compute the fixed shape
+  try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
+    if (!GetSolver()->ComputeFunction(aFunction)) {
+      SetErrorCode("Block driver failed to remove extra edges of the given shape");
+      return NULL;
+    }
+  }
+  catch (Standard_Failure) {
+    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
+    SetErrorCode(aFail->GetMessageString());
+    return NULL;
+  }
+
+  //Make a Python command
+  std::string doUnionFaces = (theOptimumNbFaces < 0) ? "False" : "True";
+  GEOM::TPythonDump(aFunction) << aCopy << " = geompy.RemoveExtraEdges("
+                               << theObject << ", " << doUnionFaces.data() << ")";
+
+  SetErrorCode(OK);
+  return aCopy;
+}
+
+//=============================================================================
+/*!
+ *  UnionFaces
+ */
+//=============================================================================
+Handle(GEOM_Object) GEOMImpl_IBlocksOperations::UnionFaces
+                                     (Handle(GEOM_Object) theObject)
+{
+  SetErrorCode(KO);
+
+  if (theObject.IsNull()) return NULL;
+
+  Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
+  if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be fixed
+
+  //Add a new Copy object
+  Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
+
+  //Add a function
+  Handle(GEOM_Function) aFunction =
+    aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_UNION_FACES);
+
+  //Check if the function is set correctly
+  if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
+
+  GEOMImpl_IBlockTrsf aTI (aFunction);
+  aTI.SetOriginal(aLastFunction);
+
+  //Compute the fixed shape
+  try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
+    if (!GetSolver()->ComputeFunction(aFunction)) {
+      SetErrorCode("Block driver failed to remove extra edges of the given shape");
+      return NULL;
+    }
+  }
+  catch (Standard_Failure) {
+    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
+    SetErrorCode(aFail->GetMessageString());
+    return NULL;
+  }
+
+  //Make a Python command
+  GEOM::TPythonDump(aFunction) << aCopy << " = geompy.UnionFaces("
+                               << theObject << ")";
+
+  SetErrorCode(OK);
+  return aCopy;
+}
+
+//=============================================================================
+/*!
+ *  CheckAndImprove
+ */
+//=============================================================================
+Handle(GEOM_Object) GEOMImpl_IBlocksOperations::CheckAndImprove
+                                             (Handle(GEOM_Object) theObject)
+{
+  SetErrorCode(KO);
+
+  if (theObject.IsNull()) return NULL;
+
+  Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
+  if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be fixed
+
+  //Add a new Copy object
+  Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
+
+  //Add a function
+  Handle(GEOM_Function) aFunction =
+    aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_COMPOUND_IMPROVE);
+
+  //Check if the function is set correctly
+  if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
+
+  GEOMImpl_IBlockTrsf aTI (aFunction);
+  aTI.SetOriginal(aLastFunction);
+
+  // -1 means do not unite faces on common surface (?except case of seam edge between them?)
+  //aTI.SetOptimumNbFaces(-1);
+  aTI.SetOptimumNbFaces(6);
+
+  //Compute the fixed shape
+  try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
+    if (!GetSolver()->ComputeFunction(aFunction)) {
+      SetErrorCode("Block driver failed to improve the given blocks compound");
+      return NULL;
+    }
+  }
+  catch (Standard_Failure) {
+    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
+    SetErrorCode(aFail->GetMessageString());
+    return NULL;
+  }
+
+  //Make a Python command
+  GEOM::TPythonDump(aFunction) << aCopy
+    << " = geompy.CheckAndImprove(" << theObject << ")";
+
+  SetErrorCode(OK);
+  return aCopy;
+}
+
 //=============================================================================
 /*!
  *  ExplodeCompoundOfBlocks
@@ -2059,7 +2780,7 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::ExplodeCompound
   Handle(GEOM_Function) aFunction;
 
   TopTools_MapOfShape mapShape;
-  TCollection_AsciiString anAsciiList = "[", anEntry;
+  TCollection_AsciiString anAsciiList, anEntry;
 
   // Map shapes
   TopTools_IndexedMapOfShape anIndices;
@@ -2068,6 +2789,9 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::ExplodeCompound
 
   // Explode
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
     for (; exp.More(); exp.Next()) {
       if (mapShape.Add(exp.Current())) {
@@ -2090,8 +2814,7 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::ExplodeCompound
 
           //Make a Python command
           TDF_Tool::Entry(anObj->GetEntry(), anEntry);
-          anAsciiList += anEntry;
-          anAsciiList += ",";
+          anAsciiList += anEntry + ", ";
         }
       }
     }
@@ -2107,23 +2830,15 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::ExplodeCompound
     return aBlocks;
   }
 
-  anAsciiList.Trunc(anAsciiList.Length() - 1);
-  anAsciiList += "]";
+  anAsciiList.Trunc(anAsciiList.Length() - 2);
 
   //The explode doesn't change object so no new function is required.
   aFunction = theCompound->GetLastFunction();
 
   //Make a Python command
-  TCollection_AsciiString aDescr (anAsciiList);
-  aDescr += " = IBlocksOperations.ExplodeCompoundOfBlocks(";
-  TDF_Tool::Entry(theCompound->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  aDescr += TCollection_AsciiString(theMinNbFaces) + ", ";
-  aDescr += TCollection_AsciiString(theMaxNbFaces) + ")";
-
-  TCollection_AsciiString aNewDescr = aFunction->GetDescription() + "\n";
-  aNewDescr += aDescr;
-  aFunction->SetDescription(aNewDescr);
+  GEOM::TPythonDump(aFunction, /*append=*/true)
+    << "[" << anAsciiList.ToCString() << "] = geompy.MakeBlockExplode("
+    << theCompound << ", " << theMinNbFaces << ", " << theMaxNbFaces << ")";
 
   SetErrorCode(OK);
   return aBlocks;
@@ -2153,70 +2868,83 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockNearPoint
   }
   if (aBlockOrComp.ShapeType() != TopAbs_COMPOUND &&
       aBlockOrComp.ShapeType() != TopAbs_COMPSOLID) {
-    SetErrorCode("Shape is neither a block, nor a compound of blocks");
+    SetErrorCode("Shape to find block in is not a compound");
     return NULL;
   }
 
   TopoDS_Shape anArg = thePoint->GetValue();
   if (anArg.IsNull()) {
-    SetErrorCode("Null shape is given as argument");
+    SetErrorCode("Point is null");
     return NULL;
   }
   if (anArg.ShapeType() != TopAbs_VERTEX) {
-    SetErrorCode("Element for block identification is not a vertex");
+    SetErrorCode("Shape for block identification is not a vertex");
     return NULL;
   }
 
   //Compute the Block value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     TopoDS_Shape aShape;
+
     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
     gp_Pnt aPnt = BRep_Tool::Pnt(aVert);
+    Standard_Real PX, PY, PZ;
+    aPnt.Coord(PX, PY, PZ);
 
-    // 1. Explode compound on blocks
-    TopTools_MapOfShape mapShape;
-    Standard_Integer nbSolids = 0;
+    // 1. Classify the point relatively each block
+    Standard_Integer nearest = 2, nbFound = 0;
+    TopTools_DataMapOfShapeInteger mapShapeDist;
     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
     for (; exp.More(); exp.Next()) {
-      if (mapShape.Add(exp.Current())) {
-        nbSolids++;
-      }
-    }
+      TopoDS_Shape aSolid = exp.Current();
+
+      if (!mapShapeDist.IsBound(aSolid)) {
+        Standard_Integer aDistance = 2;
+
+        // 1.a. Classify relatively Bounding box
+        Standard_Real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
+        Bnd_Box BB;
+        BRepBndLib::Add(aSolid, BB);
+        BB.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
+        if (PX < Xmin || Xmax < PX ||
+            PY < Ymin || Ymax < PY ||
+            PZ < Zmin || Zmax < PZ) {
+          // OUT of bounding box
+          aDistance = 1;
+        } else {
+          // 1.b. Classify relatively the solid itself
+          BRepClass3d_SolidClassifier SC (aSolid, aPnt, Precision::Confusion());
+          if (SC.State() == TopAbs_IN) {
+            aDistance = -1;
+          } else if (SC.State() == TopAbs_ON) {
+            aDistance = 0;
+          } else { // OUT
+            aDistance = 1;
+          }
+        }
 
-    mapShape.Clear();
-    Standard_Integer ind = 1;
-    TopTools_Array1OfShape aSolids (1, nbSolids);
-    TColStd_Array1OfInteger aDistances (1, nbSolids);
-    for (exp.Init(aBlockOrComp, TopAbs_SOLID); exp.More(); exp.Next()) {
-      if (mapShape.Add(exp.Current())) {
-        TopoDS_Shape aSolid = exp.Current();
-        aSolids(ind) = aSolid;
+        if (aDistance < nearest) {
+          nearest = aDistance;
+          aShape = aSolid;
+          nbFound = 1;
 
-        // 2. Classify the point relatively each block
-        BRepClass3d_SolidClassifier SC (aSolid, aPnt, Precision::Confusion());
-        if (SC.State() == TopAbs_IN) {
-          aDistances(ind) = -1;
-        } else if (SC.State() == TopAbs_ON) {
-          aDistances(ind) = 0;
-        } else { // OUT
-          aDistances(ind) = 1;
+          // A first found block, containing the point inside, will be returned.
+          // It is the solution, if there are no intersecting blocks in the compound.
+          if (nearest == -1) break;
+
+        } else if (aDistance == nearest) {
+          nbFound++;
+        } else {
         }
-        ind++;
-      }
-    }
 
-    // 3. Define block, containing the point or having minimum distance to it
-    Standard_Integer nearest = 2, nbFound = 0;
-    for (ind = 1; ind <= nbSolids; ind++) {
-      if (aDistances(ind) < nearest) {
-        nearest = aDistances(ind);
-        aShape = aSolids(ind);
-        nbFound = 1;
-      } else if (aDistances(ind) == nearest) {
-        nbFound++;
-      } else {
-      }
+        mapShapeDist.Bind(aSolid, aDistance);
+      } // if (!mapShapeDist.IsBound(aSolid))
     }
+
+    // 2. Define block, containing the point or having minimum distance to it
     if (nbFound > 1) {
       if (nearest == 0) {
         // The point is on boundary of some blocks and there are
@@ -2229,9 +2957,11 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockNearPoint
         // no blocks, having the point inside or on boundary.
         // We will get a nearest block
         Standard_Real minDist = RealLast();
-        for (ind = 1; ind <= nbSolids; ind++) {
-          if (aDistances(ind) == 1) {
-            BRepExtrema_DistShapeShape aDistTool (aVert, aSolids(ind));
+        TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
+        for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
+          if (mapShapeDistIter.Value() == 1) {
+            TopoDS_Shape aSolid = mapShapeDistIter.Key();
+            BRepExtrema_DistShapeShape aDistTool (aVert, aSolid);
             if (!aDistTool.IsDone()) {
               SetErrorCode("Can not find a distance from the given point to one of blocks");
               return NULL;
@@ -2239,29 +2969,31 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockNearPoint
             Standard_Real aDist = aDistTool.Value();
             if (aDist < minDist) {
               minDist = aDist;
-              aShape = aSolids(ind);
+              aShape = aSolid;
             }
           }
         }
       } else { // nearest == -1
-        // The point is inside some blocks.
-        // We will get a block with nearest center
-        Standard_Real minDist = RealLast();
-        for (ind = 1; ind <= nbSolids; ind++) {
-          if (aDistances(ind) == -1) {
-            GProp_GProps aSystem;
-            BRepGProp::VolumeProperties(aSolids(ind), aSystem);
-            gp_Pnt aCenterMass = aSystem.CentreOfMass();
-
-            Standard_Real aDist = aCenterMass.Distance(aPnt);
-            if (aDist < minDist) {
-              minDist = aDist;
-              aShape = aSolids(ind);
-            }
-          }
-        }
+//        // The point is inside some blocks.
+//        // We will get a block with nearest center
+//        Standard_Real minDist = RealLast();
+//        TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
+//        for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
+//          if (mapShapeDistIter.Value() == -1) {
+//            TopoDS_Shape aSolid = mapShapeDistIter.Key();
+//            GProp_GProps aSystem;
+//            BRepGProp::VolumeProperties(aSolid, aSystem);
+//            gp_Pnt aCenterMass = aSystem.CentreOfMass();
+//
+//            Standard_Real aDist = aCenterMass.Distance(aPnt);
+//            if (aDist < minDist) {
+//              minDist = aDist;
+//              aShape = aSolid;
+//            }
+//          }
+//        }
       }
-    }
+    } // if (nbFound > 1)
 
     if (nbFound == 0) {
       SetErrorCode("There are no blocks near the given point");
@@ -2280,21 +3012,11 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockNearPoint
     return NULL;
   }
 
-  //The GetBlockNearPoint() doesn't change object so no new function is required.
-  Handle(GEOM_Function) aFunction = theCompound->GetLastFunction();
+  Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aResult->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.GetBlockNearPoint(";
-  TDF_Tool::Entry(theCompound->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  TDF_Tool::Entry(thePoint->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  TCollection_AsciiString aNewDescr = aFunction->GetDescription() + "\n";
-  aNewDescr += aDescr;
-  aFunction->SetDescription(aNewDescr);
+  GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetBlockNearPoint("
+                               << theCompound << ", " << thePoint << ")";
 
   SetErrorCode(OK);
   return aResult;
@@ -2341,6 +3063,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockByParts
 
   //Compute the Block value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     // 1. Explode compound on solids
     TopTools_MapOfShape mapShape;
     Standard_Integer nbSolids = 0;
@@ -2404,20 +3129,11 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockByParts
     return NULL;
   }
 
-  //The GetBlockByParts() doesn't change object so no new function is required.
-  Handle(GEOM_Function) aFunction = theCompound->GetLastFunction();
+  Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
 
   //Make a Python command
-  TDF_Tool::Entry(aResult->GetEntry(), anEntry);
-  TCollection_AsciiString aDescr (anEntry);
-  aDescr += " = IBlocksOperations.GetBlockByParts(";
-  TDF_Tool::Entry(theCompound->GetEntry(), anEntry);
-  aDescr += anEntry + ", [";
-  aDescr += aPartsDescr + "])";
-
-  TCollection_AsciiString aNewDescr = aFunction->GetDescription() + "\n";
-  aNewDescr += aDescr;
-  aFunction->SetDescription(aNewDescr);
+  GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetBlockByParts("
+    << theCompound << ", [" << aPartsDescr.ToCString() << "])";
 
   SetErrorCode(OK);
   return aResult;
@@ -2445,7 +3161,8 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::GetBlocksByPart
   //Get the parts
   Standard_Integer argi, aLen = theParts->Length();
   TopTools_Array1OfShape anArgs (1, aLen);
-  TCollection_AsciiString anEntry, aPartsDescr, anAsciiList = "[";
+  TCollection_AsciiString anEntry, aPartsDescr, anAsciiList;
+
   for (argi = 1; argi <= aLen; argi++) {
     Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast(theParts->Value(argi));
     Handle(GEOM_Function) aRef = anObj->GetLastFunction();
@@ -2460,12 +3177,14 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::GetBlocksByPart
 
     // For Python command
     TDF_Tool::Entry(anObj->GetEntry(), anEntry);
-    if (argi > 1) aPartsDescr += ", ";
-    aPartsDescr += anEntry;
+    aPartsDescr += anEntry + ", ";
   }
 
   //Get the Blocks
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     TopTools_MapOfShape mapShape;
     Standard_Integer nbSolids = 0;
     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
@@ -2524,10 +3243,11 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::GetBlocksByPart
         anObj = GetEngine()->AddSubShape(theCompound, anArray);
         aBlocks->Append(anObj);
 
-        //Make a Python command
+        // For Python command
         TDF_Tool::Entry(anObj->GetEntry(), anEntry);
-        anAsciiList += anEntry;
-        anAsciiList += ",";
+        anAsciiList += anEntry + ", ";
+        if (aFunction.IsNull())
+          aFunction = anObj->GetLastFunction();
       }
     }
   }
@@ -2537,22 +3257,13 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::GetBlocksByPart
     return NULL;
   }
 
-  anAsciiList.Trunc(anAsciiList.Length() - 1);
-  anAsciiList += "]";
-
-  //The GetBlocksByParts() doesn't change object so no new function is required.
-  aFunction = theCompound->GetLastFunction();
-
   //Make a Python command
-  TCollection_AsciiString aDescr (anAsciiList);
-  aDescr += " = IBlocksOperations.GetBlocksByParts(";
-  TDF_Tool::Entry(theCompound->GetEntry(), anEntry);
-  aDescr += anEntry + ", [";
-  aDescr += aPartsDescr + "])";
+  aPartsDescr.Trunc(aPartsDescr.Length() - 2);
+  anAsciiList.Trunc(anAsciiList.Length() - 2);
 
-  TCollection_AsciiString aNewDescr = aFunction->GetDescription() + "\n";
-  aNewDescr += aDescr;
-  aFunction->SetDescription(aNewDescr);
+  GEOM::TPythonDump(aFunction) << "[" << anAsciiList.ToCString()
+    << "] = geompy.GetBlocksByParts(" << theCompound
+      << ", [" << aPartsDescr.ToCString() << "])";
 
   SetErrorCode(OK);
   return aBlocks;
@@ -2594,6 +3305,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation1D
 
   //Compute the transformation
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to make multi-transformation");
       return NULL;
@@ -2606,16 +3320,8 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation1D
   }
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aCopy->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.MakeMultiTransformation1D(";
-  TDF_Tool::Entry(theObject->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  aDescr += TCollection_AsciiString(theDirFace1) + ", ";
-  aDescr += TCollection_AsciiString(theDirFace2) + ", ";
-  aDescr += TCollection_AsciiString(theNbTimes)  + ") ";
-
-  aFunction->SetDescription(aDescr);
+  GEOM::TPythonDump(aFunction) << aCopy << " = geompy.MakeMultiTransformation1D("
+    << theObject << ", " << theDirFace1 << ", " << theDirFace2 << ", " << theNbTimes << ")";
 
   SetErrorCode(OK);
   return aCopy;
@@ -2663,6 +3369,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation2D
 
   //Compute the transformation
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to make multi-transformation");
       return NULL;
@@ -2675,20 +3384,165 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation2D
   }
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aCopy->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.MakeMultiTransformation2D(";
-  TDF_Tool::Entry(theObject->GetEntry(), anEntry);
-  aDescr += anEntry + ", ";
-  aDescr += TCollection_AsciiString(theDirFace1U) + ", ";
-  aDescr += TCollection_AsciiString(theDirFace2U) + ", ";
-  aDescr += TCollection_AsciiString(theNbTimesU)  + ", ";
-  aDescr += TCollection_AsciiString(theDirFace1V) + ", ";
-  aDescr += TCollection_AsciiString(theDirFace2V) + ", ";
-  aDescr += TCollection_AsciiString(theNbTimesV)  + ") ";
-
-  aFunction->SetDescription(aDescr);
+  GEOM::TPythonDump(aFunction) << aCopy << " = geompy.MakeMultiTransformation2D("
+    << theObject << ", " << theDirFace1U << ", " << theDirFace2U << ", " << theNbTimesU
+      << ", " << theDirFace1V << ", " << theDirFace2V << ", " << theNbTimesV << ")";
 
   SetErrorCode(OK);
   return aCopy;
 }
+
+//=============================================================================
+/*!
+ *  Propagate
+ */
+//=============================================================================
+Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::Propagate
+                                                 (Handle(GEOM_Object) theShape)
+{
+  SetErrorCode(KO);
+
+  if (theShape.IsNull()) return NULL;
+
+  TopoDS_Shape aShape = theShape->GetValue();
+  if (aShape.IsNull()) return NULL;
+
+  TopTools_IndexedMapOfShape anIndices;
+  TopExp::MapShapes(aShape, anIndices);
+
+  TopTools_IndexedDataMapOfShapeListOfShape MEW;
+  GEOMImpl_Block6Explorer::MapShapesAndAncestors
+    (aShape, TopAbs_EDGE, TopAbs_WIRE, MEW);
+  Standard_Integer ie, nbEdges = MEW.Extent();
+
+  // Result
+  Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
+
+  TopTools_MapOfShape mapAcceptedEdges;
+  TCollection_AsciiString aListRes, anEntry;
+
+  // Sort shapes in current chain (Mantis issue 21053)
+  TopTools_DataMapOfShapeListOfShape aMapChains;
+  TopTools_ListOfShape aFirstInChains;
+
+  for (ie = 1; ie <= nbEdges; ie++) {
+    TopoDS_Shape curE = MEW.FindKey(ie);
+
+    if (mapAcceptedEdges.Contains(curE)) continue;
+
+    // Build the chain
+    TopTools_ListOfShape currentChain;
+    TopTools_ListOfShape listPrevEdges;
+
+    currentChain.Append(curE);
+    listPrevEdges.Append(curE);
+    mapAcceptedEdges.Add(curE);
+
+    // Collect all edges pass by pass
+    while (listPrevEdges.Extent() > 0) {
+      // List of edges, added to chain on this cycle pass
+      TopTools_ListOfShape listCurEdges;
+
+      // Find the next portion of edges
+      TopTools_ListIteratorOfListOfShape itE (listPrevEdges);
+      for (; itE.More(); itE.Next()) {
+        TopoDS_Shape anE = itE.Value();
+
+        // Iterate on faces, having edge <anE>
+        TopTools_ListIteratorOfListOfShape itW (MEW.FindFromKey(anE));
+        for (; itW.More(); itW.Next()) {
+          TopoDS_Shape aW = itW.Value();
+          TopoDS_Shape anOppE;
+
+          BRepTools_WireExplorer aWE (TopoDS::Wire(aW));
+          Standard_Integer nb = 1, found = 0;
+          TopTools_Array1OfShape anEdges (1,4);
+          for (; aWE.More(); aWE.Next(), nb++) {
+            if (nb > 4) {
+              found = 0;
+              break;
+            }
+            anEdges(nb) = aWE.Current();
+            if (anEdges(nb).IsSame(anE)) found = nb;
+          }
+
+          if (nb == 5 && found > 0) {
+            // Quadrangle face found, get an opposite edge
+            Standard_Integer opp = found + 2;
+            if (opp > 4) opp -= 4;
+            anOppE = anEdges(opp);
+
+            if (!mapAcceptedEdges.Contains(anOppE)) {
+              // Add found edge to the chain
+              currentChain.Append(anOppE);
+              listCurEdges.Append(anOppE);
+              mapAcceptedEdges.Add(anOppE);
+            }
+          } // if (nb == 5 && found > 0)
+        } // for (; itF.More(); itF.Next())
+      } // for (; itE.More(); itE.Next())
+
+      listPrevEdges = listCurEdges;
+    } // while (listPrevEdges.Extent() > 0)
+
+    // Sort shapes in current chain (Mantis issue 21053)
+    GEOMUtils::SortShapes(currentChain, Standard_False);
+    aFirstInChains.Append(currentChain.First());
+    aMapChains.Bind(currentChain.First(), currentChain);
+  }
+
+  // Sort chains (Mantis issue 21053)
+  GEOMUtils::SortShapes(aFirstInChains, Standard_False);
+
+  // Store sorted chains in the document
+  TopTools_ListIteratorOfListOfShape aChainsIt (aFirstInChains);
+  for (; aChainsIt.More(); aChainsIt.Next()) {
+    TopoDS_Shape aFirstInChain = aChainsIt.Value();
+    const TopTools_ListOfShape& currentChain = aMapChains.Find(aFirstInChain);
+
+    // Store the chain in the document
+    Handle(TColStd_HArray1OfInteger) anArray =
+      new TColStd_HArray1OfInteger (1, currentChain.Extent());
+
+    // Fill array of sub-shape indices
+    TopTools_ListIteratorOfListOfShape itSub (currentChain);
+    for (int index = 1; itSub.More(); itSub.Next(), ++index) {
+      int id = anIndices.FindIndex(itSub.Value());
+      anArray->SetValue(index, id);
+    }
+
+    // Add a new group object
+    Handle(GEOM_Object) aChain = GetEngine()->AddSubShape(theShape, anArray);
+
+    // Set a GROUP type
+    aChain->SetType(GEOM_GROUP);
+
+    // Set a sub-shape type
+    TDF_Label aFreeLabel = aChain->GetFreeLabel();
+    TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)TopAbs_EDGE);
+
+    // Add the chain to the result
+    aSeq->Append(aChain);
+
+    //Make a Python command
+    TDF_Tool::Entry(aChain->GetEntry(), anEntry);
+    aListRes += anEntry + ", ";
+  }
+
+  if (aSeq->IsEmpty()) {
+    SetErrorCode("There are no quadrangle faces in the shape");
+    return aSeq;
+  }
+
+  aListRes.Trunc(aListRes.Length() - 2);
+
+  // The Propagation doesn't change object so no new function is required.
+  Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
+
+  // Make a Python command
+  GEOM::TPythonDump(aFunction, /*append=*/true)
+    << "[" << aListRes.ToCString() << "] = geompy.Propagate(" << theShape << ")";
+
+  SetErrorCode(OK);
+  return aSeq;
+}