Salome HOME
Added 2 new operations:
[modules/geom.git] / src / GEOMImpl / GEOMImpl_IBlocksOperations.cxx
index 1f433c62215c7b1cd2aeec1c1dcd40ff93bda01b..899afc083eb8e2944ef7e6d0837017870a4471b6 100644 (file)
@@ -1,26 +1,55 @@
-using namespace std;
+// Copyright (C) 2007-2011  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 WNT
+#pragma warning( disable:4786 )
+#endif
+
+#include <Standard_Stream.hxx>
 
-#include "GEOMImpl_IBlocksOperations.hxx"
+#include <GEOMImpl_IBlocksOperations.hxx>
 
-#include "GEOMImpl_Types.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 <GEOMImpl_BlockDriver.hxx>
+#include <GEOMImpl_IBlocks.hxx>
+#include <GEOMImpl_IBlockTrsf.hxx>
+#include <GEOMImpl_CopyDriver.hxx>
+#include <GEOMImpl_Block6Explorer.hxx>
+#include <GEOMImpl_IShapesOperations.hxx>
 
-#include "GEOM_Function.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 <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>
@@ -70,6 +99,7 @@ using namespace std;
 
 #include <Precision.hxx>
 
+#include <Standard_Failure.hxx>
 #include <Standard_ErrorHandler.hxx> // CAREFUL ! position of this file is critic : see Lucien PIGNOLONI / OCC
 
 //=============================================================================
@@ -137,6 +167,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;
@@ -149,19 +182,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;
@@ -203,6 +225,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;
@@ -215,15 +240,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;
@@ -272,6 +290,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;
@@ -284,19 +305,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;
@@ -352,6 +362,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;
@@ -364,23 +377,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;
@@ -422,6 +421,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;
@@ -434,15 +436,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;
@@ -482,6 +477,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;
@@ -494,13 +492,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;
@@ -508,7 +501,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeBlockCompound
 
 //=============================================================================
 /*!
- *  GetEdge
+ *  GetPoint
  */
 //=============================================================================
 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetPoint
@@ -528,13 +521,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;
   }
 
@@ -576,19 +563,84 @@ 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);
+
+  // The GetPoint() 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)
+    << aResult << " = geompy.GetVertexNearPoint("
+    << theShape << ", " << thePoint << ")";
 
   SetErrorCode(OK);
   return aResult;
@@ -614,13 +666,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;
   }
 
@@ -638,6 +684,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);
@@ -696,23 +745,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;
@@ -737,13 +774,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;
   }
 
@@ -759,6 +790,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetEdgeNearPoint
 
   //Compute the Edge value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     TopoDS_Shape aShape;
 
     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
@@ -773,6 +807,11 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetEdgeNearPoint
       }
     }
 
+    if (nbEdges == 0) {
+      SetErrorCode("Given shape contains no edges");
+      return NULL;
+    }
+
     mapShape.Clear();
     Standard_Integer ind = 1;
     TopTools_Array1OfShape anEdges (1, nbEdges);
@@ -826,21 +865,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;
@@ -873,12 +902,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();
@@ -899,6 +922,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;
@@ -982,27 +1008,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;
@@ -1031,12 +1042,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();
@@ -1052,6 +1057,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;
@@ -1125,23 +1133,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;
@@ -1186,6 +1182,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;
@@ -1204,21 +1203,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;
@@ -1233,9 +1222,6 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceNearPoint
                                                 (Handle(GEOM_Object) theShape,
                                                  Handle(GEOM_Object) thePoint)
 {
-//  OSD_Timer timer1, timer2, timer3, timer4, timer5;
-//  timer1.Start();
-
   SetErrorCode(KO);
 
   //New object
@@ -1249,12 +1235,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()) {
@@ -1268,6 +1248,9 @@ 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);
@@ -1275,9 +1258,6 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceNearPoint
     Standard_Real PX, PY, PZ;
     aPnt.Coord(PX, PY, PZ);
 
-//    timer1.Stop();
-//    timer2.Start();
-
     // 1. Classify the point relatively each face
     Standard_Integer nearest = 2, nbFound = 0;
     TopTools_DataMapOfShapeInteger mapShapeDist;
@@ -1328,9 +1308,6 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceNearPoint
       } // if (!mapShapeDist.IsBound(aFace))
     }
 
-//    timer2.Stop();
-//    timer3.Start();
-
     // 2. Define face, containing the point or having minimum distance to it
     if (nbFound > 1) {
       if (nearest == 0) {
@@ -1401,9 +1378,6 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceNearPoint
       }
     } // if (nbFound > 1)
 
-//    timer3.Stop();
-//    timer4.Start();
-
     if (nbFound == 0) {
       SetErrorCode("There are no faces near the given point");
       return NULL;
@@ -1414,8 +1388,6 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceNearPoint
       anArray->SetValue(1, anIndices.FindIndex(aShape));
       aResult = GetEngine()->AddSubShape(theShape, anArray);
     }
-
-//    timer4.Stop();
   }
   catch (Standard_Failure) {
     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
@@ -1423,35 +1395,13 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceNearPoint
     return NULL;
   }
 
-//  timer5.Start();
-
-  //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);
-
-//  timer5.Stop();
-//
-//  cout << "Show current face times:" << endl;
-//  timer1.Show();
-//  timer2.Show();
-//  timer3.Show();
-//  timer4.Show();
-//  timer5.Show();
-
   return aResult;
 }
 
@@ -1477,12 +1427,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()) {
@@ -1496,6 +1440,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);
@@ -1578,21 +1525,142 @@ 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
+  GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceByNormale("
+    << theShape << ", " << theVector << ")";
+
+  SetErrorCode(OK);
+  return aResult;
+}
+
+//=============================================================================
+/*!
+ *  GetShapesNearPoint
+ */
+//=============================================================================
+Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetShapesNearPoint
+                                         (Handle(GEOM_Object)    theShape,
+                                          Handle(GEOM_Object)    thePoint,
+                                          const Standard_Integer theShapeType,
+                                          const Standard_Real    theConstTolerance)
+{
+  SetErrorCode(KO);
+
+  // New object
+  Handle(GEOM_Object) aResult;
+
+  // 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;
+    Standard_Integer nbEdges = 0;
+    TopExp_Explorer exp (aBlockOrComp, TopAbs_ShapeEnum(theShapeType));
+    for (; exp.More(); exp.Next()) {
+      if (mapShape.Add(exp.Current())) {
+        nbEdges++;
+      }
+    }
+
+    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;
+        }
+        aDistances(ind) = aDistTool.Value();
+        if (aDistances(ind) < aMinDist) {
+          aMinDist = aDistances(ind);
+        }
+        ind++;
+      }
+    }
+
+    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
-  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.GetShapesNearPoint(" << theShape << ", " << thePoint
+    << ", " << TopAbs_ShapeEnum(theShapeType) << ", " << theTolerance << ")";
 
   SetErrorCode(OK);
   return aResult;
@@ -1619,6 +1687,9 @@ Standard_Boolean GEOMImpl_IBlocksOperations::IsCompoundOfBlocks
   //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()) {
@@ -2063,7 +2134,7 @@ Standard_Boolean HasAnyConnection (const Standard_Integer         theBlockIndex,
 //=============================================================================
 Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocksOld
                                                 (Handle(GEOM_Object) theCompound,
-                                                 list<BCError>&      theErrors)
+                                                 std::list<BCError>&      theErrors)
 {
   SetErrorCode(KO);
 
@@ -2215,37 +2286,37 @@ Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocksOld
 //=============================================================================
 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 += "\nHexahedral solids with degenerated and/or seam edges: ";
+      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)
@@ -2264,7 +2335,7 @@ TCollection_AsciiString GEOMImpl_IBlocksOperations::PrintBCErrors
 //=============================================================================
 Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
                                               (Handle(GEOM_Object) theCompound,
-                                               list<BCError>&      theErrors)
+                                               std::list<BCError>& theErrors)
 {
   SetErrorCode(KO);
 
@@ -2360,7 +2431,7 @@ Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
   }
 
   // 3. Find not glued blocks
-  GEOMAlgo_GlueAnalyser aGD; 
+  GEOMAlgo_GlueAnalyser aGD;
 
   aGD.SetShape(aComp);
   aGD.SetTolerance(Precision::Confusion());
@@ -2447,7 +2518,8 @@ Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
  */
 //=============================================================================
 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::RemoveExtraEdges
-                                             (Handle(GEOM_Object) theObject)
+                                     (Handle(GEOM_Object) theObject,
+                                      const Standard_Integer theOptimumNbFaces)
 {
   SetErrorCode(KO);
 
@@ -2468,9 +2540,13 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::RemoveExtraEdges
 
   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;
@@ -2483,13 +2559,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::RemoveExtraEdges
   }
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aCopy->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.RemoveExtraEdges(";
-  TDF_Tool::Entry(theObject->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  aFunction->SetDescription(aDescr);
+  std::string doUnionFaces = (theOptimumNbFaces < 0) ? "False" : "True";
+  GEOM::TPythonDump(aFunction) << aCopy << " = geompy.RemoveExtraEdges("
+                               << theObject << ", " << doUnionFaces.data() << ")";
 
   SetErrorCode(OK);
   return aCopy;
@@ -2523,8 +2595,15 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::CheckAndImprove
   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;
@@ -2537,13 +2616,8 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::CheckAndImprove
   }
 
   //Make a Python command
-  TCollection_AsciiString anEntry, aDescr;
-  TDF_Tool::Entry(aCopy->GetEntry(), anEntry);
-  aDescr += anEntry + " = IBlocksOperations.CheckAndImprove(";
-  TDF_Tool::Entry(theObject->GetEntry(), anEntry);
-  aDescr += anEntry + ")";
-
-  aFunction->SetDescription(aDescr);
+  GEOM::TPythonDump(aFunction) << aCopy
+    << " = geompy.CheckAndImprove(" << theObject << ")";
 
   SetErrorCode(OK);
   return aCopy;
@@ -2570,7 +2644,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;
@@ -2579,6 +2653,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())) {
@@ -2601,8 +2678,7 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::ExplodeCompound
 
           //Make a Python command
           TDF_Tool::Entry(anObj->GetEntry(), anEntry);
-          anAsciiList += anEntry;
-          anAsciiList += ",";
+          anAsciiList += anEntry + ", ";
         }
       }
     }
@@ -2618,23 +2694,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;
@@ -2680,6 +2748,9 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockNearPoint
 
   //Compute the Block value
   try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
     TopoDS_Shape aShape;
 
     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
@@ -2805,21 +2876,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;
@@ -2866,6 +2927,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;
@@ -2929,20 +2993,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;
@@ -2970,7 +3025,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();
@@ -2985,12 +3041,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);
@@ -3049,10 +3107,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();
       }
     }
   }
@@ -3062,22 +3121,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;
@@ -3119,6 +3169,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;
@@ -3131,16 +3184,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;
@@ -3188,6 +3233,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;
@@ -3200,19 +3248,9 @@ 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;
@@ -3245,6 +3283,11 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::Propagate
   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);
@@ -3306,6 +3349,21 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::Propagate
       listPrevEdges = listCurEdges;
     } // while (listPrevEdges.Extent() > 0)
 
+    // Sort shapes in current chain (Mantis issue 21053)
+    GEOMImpl_IShapesOperations::SortShapes(currentChain, Standard_False);
+    aFirstInChains.Append(currentChain.First());
+    aMapChains.Bind(currentChain.First(), currentChain);
+  }
+
+  // Sort chains (Mantis issue 21053)
+  GEOMImpl_IShapesOperations::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());
@@ -3323,12 +3381,16 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::Propagate
     // Set a GROUP type
     aChain->SetType(GEOM_GROUP);
 
-    // Set a sub shape type
+    // 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()) {
@@ -3336,19 +3398,14 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::Propagate
     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
-  TCollection_AsciiString aDescr
-    ("\nlistPropagationChains = IShapesOperations.Propagate(");
-  TCollection_AsciiString anEntry;
-  TDF_Tool::Entry(theShape->GetEntry(), anEntry);
-  aDescr += (anEntry + ")");
-
-  TCollection_AsciiString anOldDescr = aFunction->GetDescription();
-  anOldDescr = anOldDescr + aDescr;
-  aFunction->SetDescription(anOldDescr);
+  GEOM::TPythonDump(aFunction, /*append=*/true)
+    << "[" << aListRes.ToCString() << "] = geompy.Propagate(" << theShape << ")";
 
   SetErrorCode(OK);
   return aSeq;