Salome HOME
Remove unnecessary includes
[modules/geom.git] / src / GEOMImpl / GEOMImpl_IBlocksOperations.cxx
index afc0223b611e73fd6061a05fee59658436738259..6a4428082015eadee67477f341e9e91d44d0d44a 100644 (file)
@@ -1,31 +1,56 @@
-using namespace std;
+// Copyright (C) 2007-2021  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, or (at your option) any later version.
+//
+// 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_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_Block6Explorer.hxx>
 
-#include "GEOM_Function.hxx"
-#include "GEOM_PythonDump.hxx"
+#include <GEOMUtils.hxx>
 
-#include "GEOMAlgo_GlueAnalyser.hxx"
-#include "GEOMAlgo_CoupleOfShapes.hxx"
-#include "GEOMAlgo_ListOfCoupleOfShapes.hxx"
-#include "GEOMAlgo_ListIteratorOfListOfCoupleOfShapes.hxx"
-#include "BlockFix_CheckTool.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 "utilities.h"
-#include "OpUtil.hxx"
-#include "Utils_ExceptHandlers.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>
 
@@ -59,6 +84,7 @@ using namespace std;
 #include <Bnd_Box.hxx>
 #include <GProp_GProps.hxx>
 
+#include <Geom_Curve.hxx>
 #include <Geom_Surface.hxx>
 #include <ShapeAnalysis_Surface.hxx>
 
@@ -71,15 +97,157 @@ 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
 
+
+/**
+ * This function returns Standard_True if the face is quadrangular. It means
+ * that it has only 1 wire with 4 edges. If there are more then 4 edges in
+ * the wire and theToleranceC1 is not negative the new implementation is used.
+ * According to it the face is quadrangular if it is quadrangular according to
+ * an old implementation or if it has a single wire with more then 4 edges
+ * that form exactly 4 bounds of C1 continuity with the given tolerance.
+ *
+ * \param theFace the face to be checked
+ * \param theToleranceC1 if negative, it is not used; otherwise it is used
+ *        to check if two neighbor edges of face have C1 continuity.
+ * \return Standard_True if the face is quadrangular; Standard_False otherwise.
+ */
+static Standard_Boolean IsQuadrangle(const TopoDS_Face   &theFace,
+                                     const Standard_Real  theToleranceC1)
+{
+  TopExp_Explorer aFExp (theFace, TopAbs_WIRE);
+
+  if (!aFExp.More()) {
+    // no wire in the face
+    return Standard_False;
+  }
+
+  TopoDS_Shape aWire = aFExp.Current();
+
+  aFExp.Next();
+
+  if (aFExp.More()) {
+    // multiple wires in the face
+    return Standard_False;
+  }
+
+  // Check number of edges in the face
+  Standard_Integer    aNbEdges = 0;
+  TopTools_MapOfShape aMapEdges;
+  TopExp_Explorer     aWExp(aWire, TopAbs_EDGE);
+
+  for (; aWExp.More(); aWExp.Next()) {
+    if (aMapEdges.Add(aWExp.Current())) {
+      aNbEdges++;
+
+      if (aNbEdges > 4) {
+        break;
+      }
+    }
+  }
+
+  if (aNbEdges < 4) {
+    return Standard_False;
+  }
+
+  if (aNbEdges > 4) {
+    if (theToleranceC1 < 0.) {
+      return Standard_False;
+    }
+
+    // Check if a wire has 4 bounds of C1 continuity.
+    BRepTools_WireExplorer aWireExp(TopoDS::Wire(aWire), theFace);
+    TopTools_ListOfShape   anEdges;
+
+    for (aNbEdges = 0; aWireExp.More(); aWireExp.Next()) {
+      const TopoDS_Edge &anEdge = aWireExp.Current();
+
+      // Skip degenerated edges.
+      if (!BRep_Tool::Degenerated(anEdge)) {
+        anEdges.Append(anEdge);
+        ++aNbEdges;
+      }
+    }
+
+    if (aNbEdges < 4) {
+      return Standard_False;
+    }
+
+    // Compute number of sharp corners.
+    anEdges.Append(anEdges.First()); // To make a loop.
+
+    TopTools_ListIteratorOfListOfShape anIter(anEdges);
+    Standard_Real      aPar[2];
+    Standard_Integer   aNbCorners = 0;
+    TopoDS_Edge        anEdge1    = TopoDS::Edge(anEdges.First());
+    Handle(Geom_Curve) aCurve1    = BRep_Tool::Curve(anEdge1, aPar[0], aPar[1]);
+    Handle(Geom_Curve) aCurve2;
+    TopoDS_Edge        anEdge2;
+    TopoDS_Vertex      aCommonVtx;
+    gp_Pnt             aPnt;
+    gp_Vec             aVec1;
+    gp_Vec             aVec2;
+    Standard_Boolean   isReversed1 = (anEdge1.Orientation() == TopAbs_REVERSED);
+    Standard_Boolean   isReversed2;
+
+    for (anIter.Next(); anIter.More(); anIter.Next()) {
+      TopoDS_Edge anEdge2 = TopoDS::Edge(anIter.Value());
+
+      if (!TopExp::CommonVertex(anEdge1, anEdge2, aCommonVtx)) {
+        // NEVERREACHED
+        return Standard_False;
+      }
+
+      // Check the angle between tangent vectors of 2 curves at this point.
+      Standard_Real aParam1 = BRep_Tool::Parameter(aCommonVtx, anEdge1);
+      Standard_Real aParam2 = BRep_Tool::Parameter(aCommonVtx, anEdge2);
+
+      aCurve2     = BRep_Tool::Curve(anEdge2, aPar[0], aPar[1]);
+      isReversed2 = (anEdge2.Orientation() == TopAbs_REVERSED);
+      aCurve1->D1(aParam1, aPnt, aVec1);
+      aCurve2->D1(aParam2, aPnt, aVec2);
+
+      if (isReversed1) {
+        aVec1.Reverse();
+      }
+
+      if (isReversed2) {
+        aVec2.Reverse();
+      }
+      const Standard_Real anAngle = aVec1.Angle(aVec2);
+
+      if (anAngle > theToleranceC1) {
+        ++aNbCorners;
+
+        if (aNbCorners > 4) {
+          break;
+        }
+      }
+
+      // Go to the next couple of edges.
+      anEdge1     = anEdge2;
+      aCurve1     = aCurve2;
+      isReversed1 = isReversed2;
+    }
+
+    // Check the total number of corners.
+    if (aNbCorners != 4) {
+      return Standard_False;
+    }
+  }
+
+  return Standard_True;
+}
+
 //=============================================================================
 /*!
  *   constructor:
  */
 //=============================================================================
-GEOMImpl_IBlocksOperations::GEOMImpl_IBlocksOperations (GEOM_Engine* theEngine, int theDocID)
-: GEOM_IOperations(theEngine, theDocID)
+GEOMImpl_IBlocksOperations::GEOMImpl_IBlocksOperations (GEOM_Engine* theEngine)
+: GEOM_IOperations(theEngine)
 {
   MESSAGE("GEOMImpl_IBlocksOperations::GEOMImpl_IBlocksOperations");
 }
@@ -110,7 +278,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad
       theEdge3.IsNull() || theEdge4.IsNull()) return NULL;
 
   //Add a new Face object
-  Handle(GEOM_Object) aFace = GetEngine()->AddObject(GetDocID(), GEOM_FACE);
+  Handle(GEOM_Object) aFace = GetEngine()->AddObject(GEOM_FACE);
 
   //Add a new Face function
   Handle(GEOM_Function) aFunction =
@@ -138,14 +306,14 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad
 
   //Compute the Face value
   try {
+    OCC_CATCH_SIGNALS;
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to compute a face");
       return NULL;
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -170,7 +338,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad2Edges
   if (theEdge1.IsNull() || theEdge2.IsNull()) return NULL;
 
   //Add a new Face object
-  Handle(GEOM_Object) aFace = GetEngine()->AddObject(GetDocID(), GEOM_FACE);
+  Handle(GEOM_Object) aFace = GetEngine()->AddObject(GEOM_FACE);
 
   //Add a new Face function
   Handle(GEOM_Function) aFunction =
@@ -193,14 +361,14 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad2Edges
 
   //Compute the Face value
   try {
+    OCC_CATCH_SIGNALS;
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to compute a face");
       return NULL;
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -227,7 +395,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad4Vertices
       thePnt3.IsNull() || thePnt4.IsNull()) return NULL;
 
   //Add a new Face object
-  Handle(GEOM_Object) aFace = GetEngine()->AddObject(GetDocID(), GEOM_FACE);
+  Handle(GEOM_Object) aFace = GetEngine()->AddObject(GEOM_FACE);
 
   //Add a new Face function
   Handle(GEOM_Function) aFunction =
@@ -255,14 +423,14 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad4Vertices
 
   //Compute the Face value
   try {
+    OCC_CATCH_SIGNALS;
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to compute a face");
       return NULL;
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -291,7 +459,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeHexa
       theFace5.IsNull() || theFace6.IsNull()) return NULL;
 
   //Add a new Solid object
-  Handle(GEOM_Object) aBlock = GetEngine()->AddObject(GetDocID(), GEOM_BLOCK);
+  Handle(GEOM_Object) aBlock = GetEngine()->AddObject(GEOM_BLOCK);
 
   //Add a new Block function
   Handle(GEOM_Function) aFunction =
@@ -324,14 +492,14 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeHexa
 
   //Compute the Block value
   try {
+    OCC_CATCH_SIGNALS;
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to compute a block");
       return NULL;
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -357,7 +525,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeHexa2Faces
   if (theFace1.IsNull() || theFace2.IsNull()) return NULL;
 
   //Add a new Solid object
-  Handle(GEOM_Object) aBlock = GetEngine()->AddObject(GetDocID(), GEOM_BLOCK);
+  Handle(GEOM_Object) aBlock = GetEngine()->AddObject(GEOM_BLOCK);
 
   //Add a new Block function
   Handle(GEOM_Function) aFunction =
@@ -380,14 +548,14 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeHexa2Faces
 
   //Compute the Block value
   try {
+    OCC_CATCH_SIGNALS;
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to compute a block");
       return NULL;
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -412,7 +580,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeBlockCompound
   if (theCompound.IsNull()) return NULL;
 
   //Add a new object
-  Handle(GEOM_Object) aBlockComp = GetEngine()->AddObject(GetDocID(), GEOM_COMPOUND);
+  Handle(GEOM_Object) aBlockComp = GetEngine()->AddObject(GEOM_COMPOUND);
 
   //Add a new BlocksComp function
   Handle(GEOM_Function) aFunction =
@@ -433,14 +601,14 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeBlockCompound
 
   //Compute the Blocks Compound value
   try {
+    OCC_CATCH_SIGNALS;
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to compute a blocks compound");
       return NULL;
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -454,7 +622,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeBlockCompound
 
 //=============================================================================
 /*!
- *  GetEdge
+ *  GetPoint
  */
 //=============================================================================
 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetPoint
@@ -474,13 +642,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;
   }
 
@@ -520,12 +682,84 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetPoint
 
   //The GetPoint() doesn't change object so no new function is required.
   Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
-  TCollection_AsciiString anOldDescr = aFunction->GetDescription();
 
   //Make a Python command
-  GEOM::TPythonDump(aFunction) << anOldDescr.ToCString() << "\n\t"
+  GEOM::TPythonDump(aFunction, /*append=*/true)
     << aResult << " = geompy.GetPoint(" << theShape << ", "
-      << theX << ", " << theY << ", " << theZ << ", " << theEpsilon << ")";
+    << 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;
@@ -551,13 +785,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;
   }
 
@@ -575,6 +803,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetEdge
 
   //Compute the Edge value
   try {
+    OCC_CATCH_SIGNALS;
     TopTools_IndexedDataMapOfShapeListOfShape MVE;
     GEOMImpl_Block6Explorer::MapShapesAndAncestors
       (aBlockOrComp, TopAbs_VERTEX, TopAbs_EDGE, MVE);
@@ -627,9 +856,8 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetEdge
       anArray->SetValue(1, anIndices.FindIndex(anEdge));
       aResult = GetEngine()->AddSubShape(theShape, anArray);
     }
-  } catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  } catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -662,13 +890,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;
   }
 
@@ -684,70 +906,18 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetEdgeNearPoint
 
   //Compute the Edge value
   try {
-    TopoDS_Shape aShape;
-
+    OCC_CATCH_SIGNALS;
     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();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -788,12 +958,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();
@@ -814,6 +978,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByPoints
 
   //Compute the Face value
   try {
+    OCC_CATCH_SIGNALS;
     TopoDS_Shape aShape;
 
     TopTools_IndexedDataMapOfShapeListOfShape MVF;
@@ -891,9 +1056,8 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByPoints
       aResult = GetEngine()->AddSubShape(theShape, anArray);
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -931,12 +1095,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();
@@ -952,6 +1110,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByEdges
 
   //Compute the Face value
   try {
+    OCC_CATCH_SIGNALS;
     TopoDS_Shape aShape;
 
     TopTools_IndexedDataMapOfShapeListOfShape MEF;
@@ -1019,9 +1178,8 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByEdges
       aResult = GetEngine()->AddSubShape(theShape, anArray);
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -1074,6 +1232,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetOppositeFace
 
   //Compute the Face value
   try {
+    OCC_CATCH_SIGNALS;
     TopoDS_Shape aShape;
 
     GEOMImpl_Block6Explorer aBlockTool;
@@ -1086,9 +1245,8 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetOppositeFace
     anArray->SetValue(1, anIndices.FindIndex(aShape));
     aResult = GetEngine()->AddSubShape(theShape, anArray);
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -1124,12 +1282,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()) {
@@ -1143,6 +1295,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceNearPoint
 
   //Compute the Face value
   try {
+    OCC_CATCH_SIGNALS;
     TopoDS_Shape aShape;
 
     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
@@ -1281,9 +1434,8 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceNearPoint
       aResult = GetEngine()->AddSubShape(theShape, anArray);
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -1319,12 +1471,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()) {
@@ -1338,6 +1484,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByNormale
 
   //Compute the Face value
   try {
+    OCC_CATCH_SIGNALS;
     TopoDS_Shape aShape;
 
     TopoDS_Edge anEdge = TopoDS::Edge(anArg);
@@ -1414,9 +1561,8 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByNormale
       aResult = GetEngine()->AddSubShape(theShape, anArray);
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -1430,6 +1576,134 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByNormale
   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 {
+    OCC_CATCH_SIGNALS;
+    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& aFail) {
+    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
@@ -1451,6 +1725,7 @@ Standard_Boolean GEOMImpl_IBlocksOperations::IsCompoundOfBlocks
   //Check
   isCompOfBlocks = Standard_True;
   try {
+    OCC_CATCH_SIGNALS;
     TopTools_MapOfShape mapShape;
     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
     for (; exp.More(); exp.Next()) {
@@ -1477,9 +1752,8 @@ Standard_Boolean GEOMImpl_IBlocksOperations::IsCompoundOfBlocks
       }
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return isCompOfBlocks;
   }
 
@@ -1495,7 +1769,9 @@ Standard_Boolean GEOMImpl_IBlocksOperations::IsCompoundOfBlocks
 void GEOMImpl_IBlocksOperations::AddBlocksFrom (const TopoDS_Shape&   theShape,
                                                 TopTools_ListOfShape& BLO,
                                                 TopTools_ListOfShape& NOT,
-                                                TopTools_ListOfShape& EXT)
+                                                TopTools_ListOfShape& EXT,
+                                                TopTools_ListOfShape& NOQ,
+                                                const Standard_Real   theToleranceC1)
 {
   TopAbs_ShapeEnum aType = theShape.ShapeType();
   switch (aType) {
@@ -1504,7 +1780,7 @@ void GEOMImpl_IBlocksOperations::AddBlocksFrom (const TopoDS_Shape&   theShape,
     {
       TopoDS_Iterator It (theShape);
       for (; It.More(); It.Next()) {
-        AddBlocksFrom(It.Value(), BLO, NOT, EXT);
+        AddBlocksFrom(It.Value(), BLO, NOT, EXT, NOQ, theToleranceC1);
       }
     }
     break;
@@ -1513,6 +1789,7 @@ void GEOMImpl_IBlocksOperations::AddBlocksFrom (const TopoDS_Shape&   theShape,
       // Check, if there are seam or degenerated edges
       BlockFix_CheckTool aTool;
       aTool.SetShape(theShape);
+      aTool.SetAngTolerance(theToleranceC1);
       aTool.Perform();
       if (aTool.NbPossibleBlocks() > 0) {
         EXT.Append(theShape);
@@ -1524,38 +1801,14 @@ void GEOMImpl_IBlocksOperations::AddBlocksFrom (const TopoDS_Shape&   theShape,
         TopExp_Explorer expF (theShape, TopAbs_FACE);
 
         for (; expF.More(); expF.Next()) {
-          if (mapFaces.Add(expF.Current())) {
-            nbFaces++;
-            if (nbFaces > 6) break;
+          TopoDS_Face aF = TopoDS::Face(expF.Current());
 
-            // get wire
-            TopoDS_Shape aF = expF.Current();
-            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;
-            }
+          if (mapFaces.Add(aF)) {
+            nbFaces++;
 
-            // 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) {
+            if (!IsQuadrangle(aF, theToleranceC1)) {
               hasNonQuadr = Standard_True;
+              NOQ.Append(aF);//0021483
             }
           }
         }
@@ -1568,97 +1821,21 @@ void GEOMImpl_IBlocksOperations::AddBlocksFrom (const TopoDS_Shape&   theShape,
       }
     }
     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) {
-  case TopAbs_COMPOUND:
-  case TopAbs_COMPSOLID:
-    {
-      TopoDS_Iterator It (theShape);
-      for (; It.More(); It.Next()) {
-        AddBlocksFromOld(It.Value(), BLO, NOT, DEG, SEA);
-      }
-    }
-    break;
-  case TopAbs_SOLID:
+  case TopAbs_SHELL: //0021483
+  case TopAbs_FACE: //0021483
     {
+      // Count edges in each face
       TopTools_MapOfShape mapFaces;
       TopExp_Explorer expF (theShape, TopAbs_FACE);
-      Standard_Integer nbFaces = 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 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;
-          }
+        TopoDS_Face aF = TopoDS::Face(expF.Current());
 
-          // 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 (mapFaces.Add(aF)) {
+          if (!IsQuadrangle(aF, theToleranceC1)) {
+            NOQ.Append(aF);//0021483
           }
         }
       }
-      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);
-      }
     }
     break;
   default:
@@ -1888,170 +2065,18 @@ Standard_Boolean HasAnyConnection (const Standard_Integer         theBlockIndex,
   return Standard_False;
 }
 
-//=============================================================================
-/*!
- *  CheckCompoundOfBlocksOld
- */
-//=============================================================================
-Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocksOld
-                                                (Handle(GEOM_Object) theCompound,
-                                                 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. 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
-  AddBlocksFromOld(aBlockOrComp, BLO, NOT, DEG, SEA);
-
-  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);
-  }
-
-  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;
-    SetErrorCode(OK);
-    return isCompOfBlocks;
-  }
-  if (nbBlocks == 1) {
-    SetErrorCode(OK);
-    return isCompOfBlocks;
-  }
-
-  // Convert list of blocks into array for easy and fast access
-  Standard_Integer ibl = 1;
-  TopTools_Array1OfShape aBlocks (1, nbBlocks);
-  TopTools_ListIteratorOfListOfShape BLOit (BLO);
-  for (; BLOit.More(); BLOit.Next(), ibl++) {
-    aBlocks.SetValue(ibl, BLOit.Value());
-  }
-
-  // 2. Find relations between all blocks,
-  //    report connection errors (NOT_GLUED and INVALID_CONNECTION)
-  TColStd_Array2OfInteger aRelations (1, nbBlocks, 1, nbBlocks);
-  aRelations.Init(REL_NOT_CONNECTED);
-
-  Standard_Integer row = 1;
-  for (row = 1; row <= nbBlocks; row++) {
-    TopoDS_Shape aBlock = aBlocks.Value(row);
-
-    Standard_Integer col = row + 1;
-    for (; col <= nbBlocks; col++) {
-      Standard_Integer aRel = BlocksRelation(aBlock, aBlocks.Value(col));
-      if (aRel != REL_NOT_CONNECTED) {
-        aRelations.SetValue(row, col, aRel);
-        aRelations.SetValue(col, row, aRel);
-        if (aRel == REL_NOT_GLUED) {
-          // report connection error
-          isCompOfBlocks = Standard_False;
-          BCError anErr;
-          anErr.error = NOT_GLUED;
-          anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(row)));
-          anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(col)));
-          theErrors.push_back(anErr);
-        } else if (aRel == REL_COLLISION_VV ||
-                   aRel == REL_COLLISION_FF ||
-                   aRel == REL_COLLISION_EE ||
-                   aRel == REL_UNKNOWN) {
-          // report connection error
-          isCompOfBlocks = Standard_False;
-          BCError anErr;
-          anErr.error = INVALID_CONNECTION;
-          anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(row)));
-          anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(col)));
-          theErrors.push_back(anErr);
-        } else {
-        }
-      }
-    }
-  }
-
-  // 3. Find largest set of connected (good connection or not glued) blocks
-  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;
-      }
-    }
-  }
-
-  // 4. 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(aBlocks.Value(ibl)));
-      }
-    }
-  }
-  if (hasIsolated) {
-    isCompOfBlocks = Standard_False;
-    theErrors.push_back(anErr);
-  }
-
-  SetErrorCode(OK);
-  return isCompOfBlocks;
-}
-
 //=============================================================================
 /*!
  *  PrintBCErrors
  */
 //=============================================================================
 TCollection_AsciiString GEOMImpl_IBlocksOperations::PrintBCErrors
-                                                (Handle(GEOM_Object)  theCompound,
-                                                 const list<BCError>& theErrors)
+                                                 (Handle(GEOM_Object) /*theCompound*/,
+                                                 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;
@@ -2076,8 +2101,8 @@ TCollection_AsciiString GEOMImpl_IBlocksOperations::PrintBCErrors
       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)
@@ -2096,7 +2121,8 @@ TCollection_AsciiString GEOMImpl_IBlocksOperations::PrintBCErrors
 //=============================================================================
 Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
                                               (Handle(GEOM_Object) theCompound,
-                                               list<BCError>&      theErrors)
+                                               const Standard_Real theToleranceC1,
+                                               std::list<BCError>& theErrors)
 {
   SetErrorCode(KO);
 
@@ -2113,7 +2139,8 @@ Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
   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
-  AddBlocksFrom(aBlockOrComp, BLO, NOT, EXT);
+  TopTools_ListOfShape NOQ; // All non-quadrangular faces
+  AddBlocksFrom(aBlockOrComp, BLO, NOT, EXT, NOQ, theToleranceC1);
 
   // Report non-blocks
   if (NOT.Extent() > 0) {
@@ -2192,7 +2219,7 @@ Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
   }
 
   // 3. Find not glued blocks
-  GEOMAlgo_GlueAnalyser aGD; 
+  GEOMAlgo_GlueAnalyser aGD;
 
   aGD.SetShape(aComp);
   aGD.SetTolerance(Precision::Confusion());
@@ -2273,13 +2300,109 @@ Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
   return isCompOfBlocks;
 }
 
+//=============================================================================
+/*!
+ *  GetNonBlocks
+ */
+//=============================================================================
+Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetNonBlocks
+                                     (Handle(GEOM_Object)  theShape,
+                                      const Standard_Real  theToleranceC1,
+                                      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, theToleranceC1);
+
+  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)
+                                     (Handle(GEOM_Object) theObject,
+                                      const Standard_Integer theOptimumNbFaces)
 {
   SetErrorCode(KO);
 
@@ -2289,7 +2412,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::RemoveExtraEdges
   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);
+  Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GEOM_COPY);
 
   //Add a function
   Handle(GEOM_Function) aFunction =
@@ -2300,23 +2423,74 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::RemoveExtraEdges
 
   GEOMImpl_IBlockTrsf aTI (aFunction);
   aTI.SetOriginal(aLastFunction);
+  aTI.SetOptimumNbFaces(theOptimumNbFaces);
 
   //Compute the fixed shape
   try {
+    OCC_CATCH_SIGNALS;
     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());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
   //Make a Python command
-  GEOM::TPythonDump(aFunction) << aCopy
-    << " = geompy.RemoveExtraEdges(" << theObject << ")";
+  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(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 {
+    OCC_CATCH_SIGNALS;
+    if (!GetSolver()->ComputeFunction(aFunction)) {
+      SetErrorCode("Block driver failed to remove extra edges of the given shape");
+      return NULL;
+    }
+  }
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
+    return NULL;
+  }
+
+  //Make a Python command
+  GEOM::TPythonDump(aFunction) << aCopy << " = geompy.UnionFaces("
+                               << theObject << ")";
 
   SetErrorCode(OK);
   return aCopy;
@@ -2338,7 +2512,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::CheckAndImprove
   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);
+  Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GEOM_COPY);
 
   //Add a function
   Handle(GEOM_Function) aFunction =
@@ -2350,16 +2524,20 @@ 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 {
+    OCC_CATCH_SIGNALS;
     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());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -2401,6 +2579,7 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::ExplodeCompound
 
   // Explode
   try {
+    OCC_CATCH_SIGNALS;
     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
     for (; exp.More(); exp.Next()) {
       if (mapShape.Add(exp.Current())) {
@@ -2428,9 +2607,8 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::ExplodeCompound
       }
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return aBlocks;
   }
 
@@ -2443,12 +2621,11 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::ExplodeCompound
 
   //The explode doesn't change object so no new function is required.
   aFunction = theCompound->GetLastFunction();
-  TCollection_AsciiString anOldDescr = aFunction->GetDescription();
 
   //Make a Python command
-  GEOM::TPythonDump(aFunction) << anOldDescr.ToCString() << "\n\t["
-    << anAsciiList.ToCString() << "] = geompy.MakeBlockExplode("
-      << theCompound << ", " << theMinNbFaces << ", " << theMaxNbFaces << ")";
+  GEOM::TPythonDump(aFunction, /*append=*/true)
+    << "[" << anAsciiList.ToCString() << "] = geompy.MakeBlockExplode("
+    << theCompound << ", " << theMinNbFaces << ", " << theMaxNbFaces << ")";
 
   SetErrorCode(OK);
   return aBlocks;
@@ -2494,6 +2671,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockNearPoint
 
   //Compute the Block value
   try {
+    OCC_CATCH_SIGNALS;
     TopoDS_Shape aShape;
 
     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
@@ -2613,9 +2791,8 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockNearPoint
       aResult = GetEngine()->AddSubShape(theCompound, anArray);
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -2670,6 +2847,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockByParts
 
   //Compute the Block value
   try {
+    OCC_CATCH_SIGNALS;
     // 1. Explode compound on solids
     TopTools_MapOfShape mapShape;
     Standard_Integer nbSolids = 0;
@@ -2727,9 +2905,8 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockByParts
       anArray->SetValue(1, anIndices.FindIndex(aShape));
       aResult = GetEngine()->AddSubShape(theCompound, anArray);
     }
-  } catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  } catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -2786,6 +2963,7 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::GetBlocksByPart
 
   //Get the Blocks
   try {
+    OCC_CATCH_SIGNALS;
     TopTools_MapOfShape mapShape;
     Standard_Integer nbSolids = 0;
     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
@@ -2852,9 +3030,8 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::GetBlocksByPart
       }
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -2889,7 +3066,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation1D
   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be moved
 
   //Add a new Copy object
-  Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
+  Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GEOM_COPY);
 
   //Add a translate function
   Handle(GEOM_Function) aFunction =
@@ -2906,14 +3083,14 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation1D
 
   //Compute the transformation
   try {
+    OCC_CATCH_SIGNALS;
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to make multi-transformation");
       return NULL;
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -2947,7 +3124,7 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation2D
   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be moved
 
   //Add a new Copy object
-  Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
+  Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GEOM_COPY);
 
   //Add a translate function
   Handle(GEOM_Function) aFunction =
@@ -2967,14 +3144,14 @@ Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation2D
 
   //Compute the transformation
   try {
+    OCC_CATCH_SIGNALS;
     if (!GetSolver()->ComputeFunction(aFunction)) {
       SetErrorCode("Block driver failed to make multi-transformation");
       return NULL;
     }
   }
-  catch (Standard_Failure) {
-    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
-    SetErrorCode(aFail->GetMessageString());
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
     return NULL;
   }
 
@@ -3016,6 +3193,10 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::Propagate
   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);
 
@@ -3076,6 +3257,21 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::Propagate
       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());
@@ -3093,7 +3289,7 @@ 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);
 
@@ -3114,11 +3310,10 @@ Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::Propagate
 
   // The Propagation doesn't change object so no new function is required.
   Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
-  TCollection_AsciiString anOldDescr = aFunction->GetDescription();
 
   // Make a Python command
-  GEOM::TPythonDump(aFunction) << anOldDescr.ToCString() << "\n\t["
-    << aListRes.ToCString() << "] = geompy.Propagate(" << theShape << ")";
+  GEOM::TPythonDump(aFunction, /*append=*/true)
+    << "[" << aListRes.ToCString() << "] = geompy.Propagate(" << theShape << ")";
 
   SetErrorCode(OK);
   return aSeq;