src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

   1 // Copyright (C) 2007-2014  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #ifdef WIN32
  24 #pragma warning( disable:4786 )
  25 #endif
  26
  27 #include <Standard_Stream.hxx>
  28
  29 #include <GEOMImpl_IBlocksOperations.hxx>
  30
  31 #include <GEOMImpl_Types.hxx>
  32
  33 #include <GEOMImpl_BlockDriver.hxx>
  34 #include <GEOMImpl_IBlocks.hxx>
  35 #include <GEOMImpl_IBlockTrsf.hxx>
  36 #include <GEOMImpl_Block6Explorer.hxx>
  37
  38 #include <GEOMUtils.hxx>
  39
  40 #include <GEOM_Function.hxx>
  41 #include <GEOM_PythonDump.hxx>
  42
  43 #include <GEOMAlgo_GlueAnalyser.hxx>
  44 #include <GEOMAlgo_CoupleOfShapes.hxx>
  45 #include <GEOMAlgo_ListOfCoupleOfShapes.hxx>
  46 #include <GEOMAlgo_ListIteratorOfListOfCoupleOfShapes.hxx>
  47 #include <BlockFix_CheckTool.hxx>
  48
  49 #include <Basics_OCCTVersion.hxx>
  50
  51 #include "utilities.h"
  52 #include <OpUtil.hxx>
  53 #include <Utils_ExceptHandlers.hxx>
  54
  55 #include <TFunction_DriverTable.hxx>
  56 #include <TFunction_Driver.hxx>
  57 #include <TFunction_Logbook.hxx>
  58 #include <TDataStd_Integer.hxx>
  59 #include <TDF_Tool.hxx>
  60
  61 #include <BRep_Tool.hxx>
  62 #include <BRep_Builder.hxx>
  63 #include <BRepTools.hxx>
  64 #include <BRepTools_WireExplorer.hxx>
  65 #include <BRepGProp.hxx>
  66 #include <BRepBndLib.hxx>
  67 #include <BRepAdaptor_Surface.hxx>
  68 #include <BRepClass_FaceClassifier.hxx>
  69 #include <BRepClass3d_SolidClassifier.hxx>
  70 #include <BRepExtrema_DistShapeShape.hxx>
  71
  72 #include <TopAbs.hxx>
  73 #include <TopoDS.hxx>
  74 #include <TopoDS_Edge.hxx>
  75 #include <TopoDS_Vertex.hxx>
  76 #include <TopoDS_Compound.hxx>
  77 #include <TopoDS_Iterator.hxx>
  78 #include <TopExp.hxx>
  79 #include <TopExp_Explorer.hxx>
  80 #include <TopTools_MapOfShape.hxx>
  81 #include <TopTools_Array1OfShape.hxx>
  82 #include <TopTools_IndexedMapOfShape.hxx>
  83 #include <TopTools_DataMapOfShapeInteger.hxx>
  84 #include <TopTools_ListIteratorOfListOfShape.hxx>
  85 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
  86 #include <TopTools_DataMapIteratorOfDataMapOfShapeInteger.hxx>
  87
  88 #include <Bnd_Box.hxx>
  89 #include <GProp_GProps.hxx>
  90
  91 #include <Geom_Surface.hxx>
  92 #include <ShapeAnalysis_Surface.hxx>
  93
  94 #include <TColStd_MapOfInteger.hxx>
  95 #include <TColStd_Array1OfReal.hxx>
  96 #include <TColStd_Array1OfInteger.hxx>
  97 #include <TColStd_Array2OfInteger.hxx>
  98
  99 //#include <OSD_Timer.hxx>
 100
 101 #include <Precision.hxx>
 102
 103 #include <Standard_Failure.hxx>
 104 #include <Standard_ErrorHandler.hxx> // CAREFUL ! position of this file is critic : see Lucien PIGNOLONI / OCC
 105
 106 //=============================================================================
 107 /*!
 108  *   constructor:
 109  */
 110 //=============================================================================
 111 GEOMImpl_IBlocksOperations::GEOMImpl_IBlocksOperations (GEOM_Engine* theEngine, int theDocID)
 112 : GEOM_IOperations(theEngine, theDocID)
 113 {
 114   MESSAGE("GEOMImpl_IBlocksOperations::GEOMImpl_IBlocksOperations");
 115 }
 116
 117 //=============================================================================
 118 /*!
 119  *  destructor
 120  */
 121 //=============================================================================
 122 GEOMImpl_IBlocksOperations::~GEOMImpl_IBlocksOperations()
 123 {
 124   MESSAGE("GEOMImpl_IBlocksOperations::~GEOMImpl_IBlocksOperations");
 125 }
 126
 127
 128 //=============================================================================
 129 /*!
 130  *  MakeQuad
 131  */
 132 //=============================================================================
 133 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad
 134                      (Handle(GEOM_Object) theEdge1, Handle(GEOM_Object) theEdge2,
 135                       Handle(GEOM_Object) theEdge3, Handle(GEOM_Object) theEdge4)
 136 {
 137   SetErrorCode(KO);
 138
 139   if (theEdge1.IsNull() || theEdge2.IsNull() ||
 140       theEdge3.IsNull() || theEdge4.IsNull()) return NULL;
 141
 142   //Add a new Face object
 143   Handle(GEOM_Object) aFace = GetEngine()->AddObject(GetDocID(), GEOM_FACE);
 144
 145   //Add a new Face function
 146   Handle(GEOM_Function) aFunction =
 147     aFace->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_FACE_FOUR_EDGES);
 148
 149   //Check if the function is set correctly
 150   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
 151
 152   GEOMImpl_IBlocks aPI (aFunction);
 153
 154   Handle(GEOM_Function) aRef1 = theEdge1->GetLastFunction();
 155   Handle(GEOM_Function) aRef2 = theEdge2->GetLastFunction();
 156   Handle(GEOM_Function) aRef3 = theEdge3->GetLastFunction();
 157   Handle(GEOM_Function) aRef4 = theEdge4->GetLastFunction();
 158   if (aRef1.IsNull() || aRef2.IsNull() ||
 159       aRef3.IsNull() || aRef4.IsNull()) return NULL;
 160
 161   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
 162   aShapesSeq->Append(aRef1);
 163   aShapesSeq->Append(aRef2);
 164   aShapesSeq->Append(aRef3);
 165   aShapesSeq->Append(aRef4);
 166
 167   aPI.SetShapes(aShapesSeq);
 168
 169   //Compute the Face value
 170   try {
 171 #if OCC_VERSION_LARGE > 0x06010000
 172     OCC_CATCH_SIGNALS;
 173 #endif
 174     if (!GetSolver()->ComputeFunction(aFunction)) {
 175       SetErrorCode("Block driver failed to compute a face");
 176       return NULL;
 177     }
 178   }
 179   catch (Standard_Failure) {
 180     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 181     SetErrorCode(aFail->GetMessageString());
 182     return NULL;
 183   }
 184
 185   //Make a Python command
 186   GEOM::TPythonDump(aFunction) << aFace << " = geompy.MakeQuad("
 187     << theEdge1 << ", " << theEdge2 << ", " << theEdge3 << ", " << theEdge4 << ")";
 188
 189   SetErrorCode(OK);
 190   return aFace;
 191 }
 192
 193 //=============================================================================
 194 /*!
 195  *  MakeQuad2Edges
 196  */
 197 //=============================================================================
 198 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad2Edges
 199                      (Handle(GEOM_Object) theEdge1, Handle(GEOM_Object) theEdge2)
 200 {
 201   SetErrorCode(KO);
 202
 203   if (theEdge1.IsNull() || theEdge2.IsNull()) return NULL;
 204
 205   //Add a new Face object
 206   Handle(GEOM_Object) aFace = GetEngine()->AddObject(GetDocID(), GEOM_FACE);
 207
 208   //Add a new Face function
 209   Handle(GEOM_Function) aFunction =
 210     aFace->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_FACE_TWO_EDGES);
 211
 212   //Check if the function is set correctly
 213   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
 214
 215   GEOMImpl_IBlocks aPI (aFunction);
 216
 217   Handle(GEOM_Function) aRef1 = theEdge1->GetLastFunction();
 218   Handle(GEOM_Function) aRef2 = theEdge2->GetLastFunction();
 219   if (aRef1.IsNull() || aRef2.IsNull()) return NULL;
 220
 221   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
 222   aShapesSeq->Append(aRef1);
 223   aShapesSeq->Append(aRef2);
 224
 225   aPI.SetShapes(aShapesSeq);
 226
 227   //Compute the Face value
 228   try {
 229 #if OCC_VERSION_LARGE > 0x06010000
 230     OCC_CATCH_SIGNALS;
 231 #endif
 232     if (!GetSolver()->ComputeFunction(aFunction)) {
 233       SetErrorCode("Block driver failed to compute a face");
 234       return NULL;
 235     }
 236   }
 237   catch (Standard_Failure) {
 238     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 239     SetErrorCode(aFail->GetMessageString());
 240     return NULL;
 241   }
 242
 243   //Make a Python command
 244   GEOM::TPythonDump(aFunction) << aFace << " = geompy.MakeQuad2Edges("
 245                                << theEdge1 << ", " << theEdge2 << ")";
 246
 247   SetErrorCode(OK);
 248   return aFace;
 249 }
 250
 251 //=============================================================================
 252 /*!
 253  *  MakeQuad4Vertices
 254  */
 255 //=============================================================================
 256 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad4Vertices
 257                      (Handle(GEOM_Object) thePnt1, Handle(GEOM_Object) thePnt2,
 258                       Handle(GEOM_Object) thePnt3, Handle(GEOM_Object) thePnt4)
 259 {
 260   SetErrorCode(KO);
 261
 262   if (thePnt1.IsNull() || thePnt2.IsNull() ||
 263       thePnt3.IsNull() || thePnt4.IsNull()) return NULL;
 264
 265   //Add a new Face object
 266   Handle(GEOM_Object) aFace = GetEngine()->AddObject(GetDocID(), GEOM_FACE);
 267
 268   //Add a new Face function
 269   Handle(GEOM_Function) aFunction =
 270     aFace->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_FACE_FOUR_PNT);
 271
 272   //Check if the function is set correctly
 273   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
 274
 275   GEOMImpl_IBlocks aPI (aFunction);
 276
 277   Handle(GEOM_Function) aRef1 = thePnt1->GetLastFunction();
 278   Handle(GEOM_Function) aRef2 = thePnt2->GetLastFunction();
 279   Handle(GEOM_Function) aRef3 = thePnt3->GetLastFunction();
 280   Handle(GEOM_Function) aRef4 = thePnt4->GetLastFunction();
 281   if (aRef1.IsNull() || aRef2.IsNull() ||
 282       aRef3.IsNull() || aRef4.IsNull()) return NULL;
 283
 284   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
 285   aShapesSeq->Append(aRef1);
 286   aShapesSeq->Append(aRef2);
 287   aShapesSeq->Append(aRef3);
 288   aShapesSeq->Append(aRef4);
 289
 290   aPI.SetShapes(aShapesSeq);
 291
 292   //Compute the Face value
 293   try {
 294 #if OCC_VERSION_LARGE > 0x06010000
 295     OCC_CATCH_SIGNALS;
 296 #endif
 297     if (!GetSolver()->ComputeFunction(aFunction)) {
 298       SetErrorCode("Block driver failed to compute a face");
 299       return NULL;
 300     }
 301   }
 302   catch (Standard_Failure) {
 303     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 304     SetErrorCode(aFail->GetMessageString());
 305     return NULL;
 306   }
 307
 308   //Make a Python command
 309   GEOM::TPythonDump(aFunction) << aFace << " = geompy.MakeQuad4Vertices("
 310     << thePnt1 << ", " << thePnt2 << ", " << thePnt3 << ", " << thePnt4 << ")";
 311
 312   SetErrorCode(OK);
 313   return aFace;
 314 }
 315
 316 //=============================================================================
 317 /*!
 318  *  MakeHexa
 319  */
 320 //=============================================================================
 321 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeHexa
 322                      (Handle(GEOM_Object) theFace1, Handle(GEOM_Object) theFace2,
 323                       Handle(GEOM_Object) theFace3, Handle(GEOM_Object) theFace4,
 324                       Handle(GEOM_Object) theFace5, Handle(GEOM_Object) theFace6)
 325 {
 326   SetErrorCode(KO);
 327
 328   if (theFace1.IsNull() || theFace2.IsNull() ||
 329       theFace3.IsNull() || theFace4.IsNull() ||
 330       theFace5.IsNull() || theFace6.IsNull()) return NULL;
 331
 332   //Add a new Solid object
 333   Handle(GEOM_Object) aBlock = GetEngine()->AddObject(GetDocID(), GEOM_BLOCK);
 334
 335   //Add a new Block function
 336   Handle(GEOM_Function) aFunction =
 337     aBlock->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_SIX_FACES);
 338
 339   //Check if the function is set correctly
 340   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
 341
 342   GEOMImpl_IBlocks aPI (aFunction);
 343
 344   Handle(GEOM_Function) aRef1 = theFace1->GetLastFunction();
 345   Handle(GEOM_Function) aRef2 = theFace2->GetLastFunction();
 346   Handle(GEOM_Function) aRef3 = theFace3->GetLastFunction();
 347   Handle(GEOM_Function) aRef4 = theFace4->GetLastFunction();
 348   Handle(GEOM_Function) aRef5 = theFace5->GetLastFunction();
 349   Handle(GEOM_Function) aRef6 = theFace6->GetLastFunction();
 350   if (aRef1.IsNull() || aRef2.IsNull() ||
 351       aRef3.IsNull() || aRef4.IsNull() ||
 352       aRef5.IsNull() || aRef6.IsNull()) return NULL;
 353
 354   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
 355   aShapesSeq->Append(aRef1);
 356   aShapesSeq->Append(aRef2);
 357   aShapesSeq->Append(aRef3);
 358   aShapesSeq->Append(aRef4);
 359   aShapesSeq->Append(aRef5);
 360   aShapesSeq->Append(aRef6);
 361
 362   aPI.SetShapes(aShapesSeq);
 363
 364   //Compute the Block value
 365   try {
 366 #if OCC_VERSION_LARGE > 0x06010000
 367     OCC_CATCH_SIGNALS;
 368 #endif
 369     if (!GetSolver()->ComputeFunction(aFunction)) {
 370       SetErrorCode("Block driver failed to compute a block");
 371       return NULL;
 372     }
 373   }
 374   catch (Standard_Failure) {
 375     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 376     SetErrorCode(aFail->GetMessageString());
 377     return NULL;
 378   }
 379
 380   //Make a Python command
 381   GEOM::TPythonDump(aFunction) << aBlock << " = geompy.MakeHexa("
 382     << theFace1 << ", " << theFace2 << ", " << theFace3 << ", "
 383       << theFace4 << ", " << theFace5 << ", " << theFace6 << ")";
 384
 385   SetErrorCode(OK);
 386   return aBlock;
 387 }
 388
 389 //=============================================================================
 390 /*!
 391  *  MakeHexa2Faces
 392  */
 393 //=============================================================================
 394 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeHexa2Faces
 395                    (Handle(GEOM_Object) theFace1, Handle(GEOM_Object) theFace2)
 396 {
 397   SetErrorCode(KO);
 398
 399   if (theFace1.IsNull() || theFace2.IsNull()) return NULL;
 400
 401   //Add a new Solid object
 402   Handle(GEOM_Object) aBlock = GetEngine()->AddObject(GetDocID(), GEOM_BLOCK);
 403
 404   //Add a new Block function
 405   Handle(GEOM_Function) aFunction =
 406     aBlock->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_TWO_FACES);
 407
 408   //Check if the function is set correctly
 409   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
 410
 411   GEOMImpl_IBlocks aPI (aFunction);
 412
 413   Handle(GEOM_Function) aRef1 = theFace1->GetLastFunction();
 414   Handle(GEOM_Function) aRef2 = theFace2->GetLastFunction();
 415   if (aRef1.IsNull() || aRef2.IsNull()) return NULL;
 416
 417   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
 418   aShapesSeq->Append(aRef1);
 419   aShapesSeq->Append(aRef2);
 420
 421   aPI.SetShapes(aShapesSeq);
 422
 423   //Compute the Block value
 424   try {
 425 #if OCC_VERSION_LARGE > 0x06010000
 426     OCC_CATCH_SIGNALS;
 427 #endif
 428     if (!GetSolver()->ComputeFunction(aFunction)) {
 429       SetErrorCode("Block driver failed to compute a block");
 430       return NULL;
 431     }
 432   }
 433   catch (Standard_Failure) {
 434     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 435     SetErrorCode(aFail->GetMessageString());
 436     return NULL;
 437   }
 438
 439   //Make a Python command
 440   GEOM::TPythonDump(aFunction) << aBlock << " = geompy.MakeHexa2Faces("
 441                                << theFace1 << ", " << theFace2 << ")";
 442
 443   SetErrorCode(OK);
 444   return aBlock;
 445 }
 446
 447 //=============================================================================
 448 /*!
 449  *  MakeBlockCompound
 450  */
 451 //=============================================================================
 452 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeBlockCompound
 453                                               (Handle(GEOM_Object) theCompound)
 454 {
 455   SetErrorCode(KO);
 456
 457   if (theCompound.IsNull()) return NULL;
 458
 459   //Add a new object
 460   Handle(GEOM_Object) aBlockComp = GetEngine()->AddObject(GetDocID(), GEOM_COMPOUND);
 461
 462   //Add a new BlocksComp function
 463   Handle(GEOM_Function) aFunction =
 464     aBlockComp->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_COMPOUND_GLUE);
 465
 466   //Check if the function is set correctly
 467   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
 468
 469   GEOMImpl_IBlocks aPI (aFunction);
 470
 471   Handle(GEOM_Function) aRef1 = theCompound->GetLastFunction();
 472   if (aRef1.IsNull()) return NULL;
 473
 474   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
 475   aShapesSeq->Append(aRef1);
 476
 477   aPI.SetShapes(aShapesSeq);
 478
 479   //Compute the Blocks Compound value
 480   try {
 481 #if OCC_VERSION_LARGE > 0x06010000
 482     OCC_CATCH_SIGNALS;
 483 #endif
 484     if (!GetSolver()->ComputeFunction(aFunction)) {
 485       SetErrorCode("Block driver failed to compute a blocks compound");
 486       return NULL;
 487     }
 488   }
 489   catch (Standard_Failure) {
 490     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 491     SetErrorCode(aFail->GetMessageString());
 492     return NULL;
 493   }
 494
 495   //Make a Python command
 496   GEOM::TPythonDump(aFunction) << aBlockComp
 497     << " = geompy.MakeBlockCompound(" << theCompound << ")";
 498
 499   SetErrorCode(OK);
 500   return aBlockComp;
 501 }
 502
 503 //=============================================================================
 504 /*!
 505  *  GetPoint
 506  */
 507 //=============================================================================
 508 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetPoint
 509                                                (Handle(GEOM_Object) theShape,
 510                                                 const Standard_Real theX,
 511                                                 const Standard_Real theY,
 512                                                 const Standard_Real theZ,
 513                                                 const Standard_Real theEpsilon)
 514 {
 515   SetErrorCode(KO);
 516
 517   //New Point object
 518   Handle(GEOM_Object) aResult;
 519
 520   // Arguments
 521   if (theShape.IsNull()) return NULL;
 522
 523   TopoDS_Shape aBlockOrComp = theShape->GetValue();
 524   if (aBlockOrComp.IsNull()) {
 525     SetErrorCode("Given shape is null");
 526     return NULL;
 527   }
 528
 529   //Compute the Vertex value
 530   gp_Pnt P (theX, theY, theZ);
 531   Standard_Real eps = Max(theEpsilon, Precision::Confusion());
 532
 533   TopoDS_Shape V;
 534   Standard_Integer isFound = 0;
 535   TopTools_MapOfShape mapShape;
 536   TopExp_Explorer exp (aBlockOrComp, TopAbs_VERTEX);
 537
 538   for (; exp.More(); exp.Next()) {
 539     if (mapShape.Add(exp.Current())) {
 540       TopoDS_Vertex aVi = TopoDS::Vertex(exp.Current());
 541       gp_Pnt aPi = BRep_Tool::Pnt(aVi);
 542       if (aPi.Distance(P) < eps) {
 543         V = aVi;
 544         isFound++;
 545       }
 546     }
 547   }
 548
 549   if (isFound == 0) {
 550     SetErrorCode("Vertex has not been found");
 551     return NULL;
 552   } else if (isFound > 1) {
 553     SetErrorCode("Multiple vertices found by the given coordinates and epsilon");
 554     return NULL;
 555   } else {
 556     TopTools_IndexedMapOfShape anIndices;
 557     TopExp::MapShapes(aBlockOrComp, anIndices);
 558     Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
 559     anArray->SetValue(1, anIndices.FindIndex(V));
 560     aResult = GetEngine()->AddSubShape(theShape, anArray);
 561   }
 562
 563   //The GetPoint() doesn't change object so no new function is required.
 564   Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
 565
 566   //Make a Python command
 567   GEOM::TPythonDump(aFunction, /*append=*/true)
 568     << aResult << " = geompy.GetPoint(" << theShape << ", "
 569     << theX << ", " << theY << ", " << theZ << ", " << theEpsilon << ")";
 570
 571   SetErrorCode(OK);
 572   return aResult;
 573 }
 574
 575 //=============================================================================
 576 /*!
 577  *  GetVertexNearPoint
 578  */
 579 //=============================================================================
 580 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetVertexNearPoint
 581                                                (Handle(GEOM_Object) theShape,
 582                                                 Handle(GEOM_Object) thePoint)
 583 {
 584   SetErrorCode(KO);
 585
 586   // New Point object
 587   Handle(GEOM_Object) aResult;
 588
 589   // Arguments
 590   if (theShape.IsNull() || thePoint.IsNull()) return NULL;
 591
 592   TopoDS_Shape aBlockOrComp = theShape->GetValue();
 593   TopoDS_Shape aPoint       = thePoint->GetValue();
 594   if (aBlockOrComp.IsNull() || aPoint.IsNull()) {
 595     SetErrorCode("Given shape is null");
 596     return NULL;
 597   }
 598
 599   if (aPoint.ShapeType() != TopAbs_VERTEX) {
 600     SetErrorCode("Element for vertex identification is not a vertex");
 601     return NULL;
 602   }
 603
 604   TopoDS_Vertex aVert = TopoDS::Vertex(aPoint);
 605   gp_Pnt aP = BRep_Tool::Pnt(aVert);
 606
 607   // Compute the Vertex value
 608   TopoDS_Shape V;
 609   bool isFound = false;
 610   Standard_Real aDist = RealLast();
 611   TopTools_MapOfShape mapShape;
 612
 613   TopExp_Explorer exp (aBlockOrComp, TopAbs_VERTEX);
 614   for (; exp.More(); exp.Next()) {
 615     if (mapShape.Add(exp.Current())) {
 616       TopoDS_Vertex aVi = TopoDS::Vertex(exp.Current());
 617       gp_Pnt aPi = BRep_Tool::Pnt(aVi);
 618       Standard_Real aDisti = aPi.Distance(aP);
 619       if (aDisti < aDist) {
 620         V = aVi;
 621         aDist = aDisti;
 622         isFound = true;
 623       }
 624     }
 625   }
 626
 627   if (!isFound) {
 628     SetErrorCode("Vertex has not been found");
 629     return NULL;
 630   }
 631
 632   TopTools_IndexedMapOfShape anIndices;
 633   TopExp::MapShapes(aBlockOrComp, anIndices);
 634   Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
 635   anArray->SetValue(1, anIndices.FindIndex(V));
 636   aResult = GetEngine()->AddSubShape(theShape, anArray);
 637
 638   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
 639
 640   // Make a Python command
 641   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetVertexNearPoint("
 642                                << theShape << ", " << thePoint << ")";
 643
 644   SetErrorCode(OK);
 645   return aResult;
 646 }
 647
 648 //=============================================================================
 649 /*!
 650  *  GetEdge
 651  */
 652 //=============================================================================
 653 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetEdge
 654                                                 (Handle(GEOM_Object) theShape,
 655                                                  Handle(GEOM_Object) thePoint1,
 656                                                  Handle(GEOM_Object) thePoint2)
 657 {
 658   SetErrorCode(KO);
 659
 660   //New Edge object
 661   Handle(GEOM_Object) aResult;
 662
 663   // Arguments
 664   if (theShape.IsNull() || thePoint1.IsNull() || thePoint2.IsNull()) return NULL;
 665
 666   TopoDS_Shape aBlockOrComp = theShape->GetValue();
 667   if (aBlockOrComp.IsNull()) {
 668     SetErrorCode("Given shape is null");
 669     return NULL;
 670   }
 671
 672   TopoDS_Shape anArg1 = thePoint1->GetValue();
 673   TopoDS_Shape anArg2 = thePoint2->GetValue();
 674   if (anArg1.IsNull() || anArg2.IsNull()) {
 675     SetErrorCode("Null shape is given as argument");
 676     return NULL;
 677   }
 678   if (anArg1.ShapeType() != TopAbs_VERTEX ||
 679       anArg2.ShapeType() != TopAbs_VERTEX) {
 680     SetErrorCode("Element for edge identification is not a vertex");
 681     return NULL;
 682   }
 683
 684   //Compute the Edge value
 685   try {
 686 #if OCC_VERSION_LARGE > 0x06010000
 687     OCC_CATCH_SIGNALS;
 688 #endif
 689     TopTools_IndexedDataMapOfShapeListOfShape MVE;
 690     GEOMImpl_Block6Explorer::MapShapesAndAncestors
 691       (aBlockOrComp, TopAbs_VERTEX, TopAbs_EDGE, MVE);
 692
 693     TopoDS_Shape V1,V2;
 694     Standard_Integer ish, ext = MVE.Extent();
 695
 696     if (MVE.Contains(anArg1)) {
 697       V1 = anArg1;
 698     } else {
 699       for (ish = 1; ish <= ext; ish++) {
 700         TopoDS_Shape aShi = MVE.FindKey(ish);
 701         if (BRepTools::Compare(TopoDS::Vertex(anArg1), TopoDS::Vertex(aShi))) {
 702           V1 = aShi;
 703           break;
 704         }
 705       }
 706     }
 707
 708     if (MVE.Contains(anArg2)) {
 709       V2 = anArg2;
 710     } else {
 711       for (ish = 1; ish <= ext; ish++) {
 712         TopoDS_Shape aShi = MVE.FindKey(ish);
 713         if (BRepTools::Compare(TopoDS::Vertex(anArg2), TopoDS::Vertex(aShi))) {
 714           V2 = aShi;
 715           break;
 716         }
 717       }
 718     }
 719
 720     if (V1.IsNull() || V2.IsNull()) {
 721       SetErrorCode("The given vertex does not belong to the shape");
 722       return NULL;
 723     }
 724
 725     TopoDS_Shape anEdge;
 726     Standard_Integer isFound =
 727       GEOMImpl_Block6Explorer::FindEdge(anEdge, V1, V2, MVE, Standard_True);
 728     if (isFound == 0) {
 729       SetErrorCode("The given vertices do not belong to one edge of the given shape");
 730       return NULL;
 731     } else if (isFound > 1) {
 732       SetErrorCode("Multiple edges found by the given vertices of the shape");
 733       return NULL;
 734     } else {
 735       TopTools_IndexedMapOfShape anIndices;
 736       TopExp::MapShapes(aBlockOrComp, anIndices);
 737       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
 738       anArray->SetValue(1, anIndices.FindIndex(anEdge));
 739       aResult = GetEngine()->AddSubShape(theShape, anArray);
 740     }
 741   } catch (Standard_Failure) {
 742     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 743     SetErrorCode(aFail->GetMessageString());
 744     return NULL;
 745   }
 746
 747   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
 748
 749   //Make a Python command
 750   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetEdge("
 751     << theShape << ", " << thePoint1 << ", " << thePoint2 << ")";
 752
 753   SetErrorCode(OK);
 754   return aResult;
 755 }
 756
 757 //=============================================================================
 758 /*!
 759  *  GetEdgeNearPoint
 760  */
 761 //=============================================================================
 762 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetEdgeNearPoint
 763                                                 (Handle(GEOM_Object) theShape,
 764                                                  Handle(GEOM_Object) thePoint)
 765 {
 766   SetErrorCode(KO);
 767
 768   //New object
 769   Handle(GEOM_Object) aResult;
 770
 771   // Arguments
 772   if (theShape.IsNull() || thePoint.IsNull()) return NULL;
 773
 774   TopoDS_Shape aBlockOrComp = theShape->GetValue();
 775   if (aBlockOrComp.IsNull()) {
 776     SetErrorCode("Given shape is null");
 777     return NULL;
 778   }
 779
 780   TopoDS_Shape anArg = thePoint->GetValue();
 781   if (anArg.IsNull()) {
 782     SetErrorCode("Null shape is given as argument");
 783     return NULL;
 784   }
 785   if (anArg.ShapeType() != TopAbs_VERTEX) {
 786     SetErrorCode("Element for edge identification is not a vertex");
 787     return NULL;
 788   }
 789
 790   //Compute the Edge value
 791   try {
 792 #if OCC_VERSION_LARGE > 0x06010000
 793     OCC_CATCH_SIGNALS;
 794 #endif
 795     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
 796     TopoDS_Shape aShape = GEOMUtils::GetEdgeNearPoint(aBlockOrComp, aVert);
 797
 798     TopTools_IndexedMapOfShape anIndices;
 799     TopExp::MapShapes(aBlockOrComp, anIndices);
 800     Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
 801     anArray->SetValue(1, anIndices.FindIndex(aShape));
 802     aResult = GetEngine()->AddSubShape(theShape, anArray);
 803   }
 804   catch (Standard_Failure) {
 805     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 806     SetErrorCode(aFail->GetMessageString());
 807     return NULL;
 808   }
 809
 810   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
 811
 812   //Make a Python command
 813   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetEdgeNearPoint("
 814                                << theShape << ", " << thePoint << ")";
 815
 816   SetErrorCode(OK);
 817   return aResult;
 818 }
 819
 820 //=============================================================================
 821 /*!
 822  *  GetFaceByPoints
 823  */
 824 //=============================================================================
 825 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByPoints
 826                                                 (Handle(GEOM_Object) theShape,
 827                                                  Handle(GEOM_Object) thePoint1,
 828                                                  Handle(GEOM_Object) thePoint2,
 829                                                  Handle(GEOM_Object) thePoint3,
 830                                                  Handle(GEOM_Object) thePoint4)
 831 {
 832   SetErrorCode(KO);
 833
 834   //New object
 835   Handle(GEOM_Object) aResult;
 836
 837   // Arguments
 838   if (theShape.IsNull() ||
 839       thePoint1.IsNull() || thePoint2.IsNull() ||
 840       thePoint3.IsNull() || thePoint4.IsNull()) return NULL;
 841
 842   TopoDS_Shape aBlockOrComp = theShape->GetValue();
 843   if (aBlockOrComp.IsNull()) {
 844     SetErrorCode("Block or compound is null");
 845     return NULL;
 846   }
 847
 848   TopoDS_Shape anArg1 = thePoint1->GetValue();
 849   TopoDS_Shape anArg2 = thePoint2->GetValue();
 850   TopoDS_Shape anArg3 = thePoint3->GetValue();
 851   TopoDS_Shape anArg4 = thePoint4->GetValue();
 852   if (anArg1.IsNull() || anArg2.IsNull() ||
 853       anArg3.IsNull() || anArg4.IsNull()) {
 854     SetErrorCode("Null shape is given as argument");
 855     return NULL;
 856   }
 857   if (anArg1.ShapeType() != TopAbs_VERTEX ||
 858       anArg2.ShapeType() != TopAbs_VERTEX ||
 859       anArg3.ShapeType() != TopAbs_VERTEX ||
 860       anArg4.ShapeType() != TopAbs_VERTEX) {
 861     SetErrorCode("Element for face identification is not a vertex");
 862     return NULL;
 863   }
 864
 865   //Compute the Face value
 866   try {
 867 #if OCC_VERSION_LARGE > 0x06010000
 868     OCC_CATCH_SIGNALS;
 869 #endif
 870     TopoDS_Shape aShape;
 871
 872     TopTools_IndexedDataMapOfShapeListOfShape MVF;
 873     GEOMImpl_Block6Explorer::MapShapesAndAncestors(aBlockOrComp, TopAbs_VERTEX, TopAbs_FACE, MVF);
 874
 875     TopoDS_Shape V1,V2,V3,V4;
 876     Standard_Integer ish, ext = MVF.Extent();
 877
 878     if (MVF.Contains(anArg1)) {
 879       V1 = anArg1;
 880     } else {
 881       for (ish = 1; ish <= ext; ish++) {
 882         TopoDS_Shape aShi = MVF.FindKey(ish);
 883         if (BRepTools::Compare(TopoDS::Vertex(anArg1), TopoDS::Vertex(aShi))) {
 884           V1 = aShi;
 885           break;
 886         }
 887       }
 888     }
 889
 890     if (MVF.Contains(anArg2)) {
 891       V2 = anArg2;
 892     } else {
 893       for (ish = 1; ish <= ext; ish++) {
 894         TopoDS_Shape aShi = MVF.FindKey(ish);
 895         if (BRepTools::Compare(TopoDS::Vertex(anArg2), TopoDS::Vertex(aShi))) {
 896           V2 = aShi;
 897           break;
 898         }
 899       }
 900     }
 901
 902     if (MVF.Contains(anArg3)) {
 903       V3 = anArg3;
 904     } else {
 905       for (ish = 1; ish <= ext; ish++) {
 906         TopoDS_Shape aShi = MVF.FindKey(ish);
 907         if (BRepTools::Compare(TopoDS::Vertex(anArg3), TopoDS::Vertex(aShi))) {
 908           V3 = aShi;
 909           break;
 910         }
 911       }
 912     }
 913
 914     if (MVF.Contains(anArg4)) {
 915       V4 = anArg4;
 916     } else {
 917       for (ish = 1; ish <= ext; ish++) {
 918         TopoDS_Shape aShi = MVF.FindKey(ish);
 919         if (BRepTools::Compare(TopoDS::Vertex(anArg4), TopoDS::Vertex(aShi))) {
 920           V4 = aShi;
 921           break;
 922         }
 923       }
 924     }
 925
 926     if (V1.IsNull() || V2.IsNull() || V3.IsNull() || V4.IsNull()) {
 927       SetErrorCode("The given vertex does not belong to the shape");
 928       return NULL;
 929     }
 930
 931     Standard_Integer isFound =
 932       GEOMImpl_Block6Explorer::FindFace(aShape, V1, V2, V3, V4, MVF, Standard_True);
 933     if (isFound == 0) {
 934       SetErrorCode("The given vertices do not belong to one face of the given shape");
 935       return NULL;
 936     } else if (isFound > 1) {
 937       SetErrorCode("The given vertices belong to several faces of the given shape");
 938       return NULL;
 939     } else {
 940       TopTools_IndexedMapOfShape anIndices;
 941       TopExp::MapShapes(aBlockOrComp, anIndices);
 942       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
 943       anArray->SetValue(1, anIndices.FindIndex(aShape));
 944       aResult = GetEngine()->AddSubShape(theShape, anArray);
 945     }
 946   }
 947   catch (Standard_Failure) {
 948     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 949     SetErrorCode(aFail->GetMessageString());
 950     return NULL;
 951   }
 952
 953   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
 954
 955   //Make a Python command
 956   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceByPoints("
 957     << theShape << ", " << thePoint1 << ", " << thePoint2
 958       << ", " << thePoint3 << ", " << thePoint4 << ")";
 959
 960   SetErrorCode(OK);
 961   return aResult;
 962 }
 963
 964 //=============================================================================
 965 /*!
 966  *  GetFaceByEdges
 967  */
 968 //=============================================================================
 969 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByEdges
 970                                                 (Handle(GEOM_Object) theShape,
 971                                                  Handle(GEOM_Object) theEdge1,
 972                                                  Handle(GEOM_Object) theEdge2)
 973 {
 974   SetErrorCode(KO);
 975
 976   //New object
 977   Handle(GEOM_Object) aResult;
 978
 979   // Arguments
 980   if (theShape.IsNull() || theEdge1.IsNull() || theEdge2.IsNull()) return NULL;
 981
 982   TopoDS_Shape aBlockOrComp = theShape->GetValue();
 983   if (aBlockOrComp.IsNull()) {
 984     SetErrorCode("Block or compound is null");
 985     return NULL;
 986   }
 987
 988   TopoDS_Shape anArg1 = theEdge1->GetValue();
 989   TopoDS_Shape anArg2 = theEdge2->GetValue();
 990   if (anArg1.IsNull() || anArg2.IsNull()) {
 991     SetErrorCode("Null shape is given as argument");
 992     return NULL;
 993   }
 994   if (anArg1.ShapeType() != TopAbs_EDGE ||
 995       anArg2.ShapeType() != TopAbs_EDGE) {
 996     SetErrorCode("Element for face identification is not an edge");
 997     return NULL;
 998   }
 999
1000   //Compute the Face value
1001   try {
1002 #if OCC_VERSION_LARGE > 0x06010000
1003     OCC_CATCH_SIGNALS;
1004 #endif
1005     TopoDS_Shape aShape;
1006
1007     TopTools_IndexedDataMapOfShapeListOfShape MEF;
1008     GEOMImpl_Block6Explorer::MapShapesAndAncestors(aBlockOrComp, TopAbs_EDGE, TopAbs_FACE, MEF);
1009
1010     TopoDS_Shape E1,E2;
1011     Standard_Integer ish, ext = MEF.Extent();
1012
1013     if (MEF.Contains(anArg1)) {
1014       E1 = anArg1;
1015     } else {
1016       for (ish = 1; ish <= ext; ish++) {
1017         TopoDS_Shape aShi = MEF.FindKey(ish);
1018         if (GEOMImpl_Block6Explorer::IsSimilarEdges(anArg1, aShi)) {
1019           E1 = aShi;
1020         }
1021       }
1022     }
1023
1024     if (MEF.Contains(anArg2)) {
1025       E2 = anArg2;
1026     } else {
1027       for (ish = 1; ish <= ext; ish++) {
1028         TopoDS_Shape aShi = MEF.FindKey(ish);
1029         if (GEOMImpl_Block6Explorer::IsSimilarEdges(anArg2, aShi)) {
1030           E2 = aShi;
1031         }
1032       }
1033     }
1034
1035     if (E1.IsNull() || E2.IsNull()) {
1036       SetErrorCode("The given edge does not belong to the shape");
1037       return NULL;
1038     }
1039
1040     const TopTools_ListOfShape& aFacesOfE1 = MEF.FindFromKey(E1);
1041     const TopTools_ListOfShape& aFacesOfE2 = MEF.FindFromKey(E2);
1042
1043     Standard_Integer isFound = 0;
1044     TopTools_ListIteratorOfListOfShape anIterF1 (aFacesOfE1);
1045     for (; anIterF1.More(); anIterF1.Next()) {
1046
1047       TopTools_ListIteratorOfListOfShape anIterF2 (aFacesOfE2);
1048       for (; anIterF2.More(); anIterF2.Next()) {
1049
1050         if (anIterF1.Value().IsSame(anIterF2.Value())) {
1051           isFound++;
1052
1053           // Store the face, defined by two edges
1054           aShape = anIterF1.Value();
1055         }
1056       }
1057     }
1058     if (isFound == 0) {
1059       SetErrorCode("The given edges do not belong to one face of the given shape");
1060       return NULL;
1061     } else if (isFound > 1) {
1062       SetErrorCode("The given edges belong to several faces of the given shape");
1063       return NULL;
1064     } else {
1065       TopTools_IndexedMapOfShape anIndices;
1066       TopExp::MapShapes(aBlockOrComp, anIndices);
1067       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
1068       anArray->SetValue(1, anIndices.FindIndex(aShape));
1069       aResult = GetEngine()->AddSubShape(theShape, anArray);
1070     }
1071   }
1072   catch (Standard_Failure) {
1073     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1074     SetErrorCode(aFail->GetMessageString());
1075     return NULL;
1076   }
1077
1078   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1079
1080   //Make a Python command
1081   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceByEdges("
1082     << theShape << ", " << theEdge1 << ", " << theEdge2 << ")";
1083
1084   SetErrorCode(OK);
1085   return aResult;
1086 }
1087
1088 //=============================================================================
1089 /*!
1090  *  GetOppositeFace
1091  */
1092 //=============================================================================
1093 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetOppositeFace
1094                                                 (Handle(GEOM_Object) theShape,
1095                                                  Handle(GEOM_Object) theFace)
1096 {
1097   SetErrorCode(KO);
1098
1099   //New object
1100   Handle(GEOM_Object) aResult;
1101
1102   // Arguments
1103   if (theShape.IsNull() || theFace.IsNull()) return NULL;
1104
1105   TopoDS_Shape aBlockOrComp = theShape->GetValue();
1106   if (aBlockOrComp.IsNull()) {
1107     SetErrorCode("Block is null");
1108     return NULL;
1109   }
1110   if (aBlockOrComp.ShapeType() != TopAbs_SOLID) {
1111     SetErrorCode("Shape is not a block");
1112     return NULL;
1113   }
1114
1115   TopoDS_Shape anArg = theFace->GetValue();
1116   if (anArg.IsNull()) {
1117     SetErrorCode("Null shape is given as argument");
1118     return NULL;
1119   }
1120   if (anArg.ShapeType() != TopAbs_FACE) {
1121     SetErrorCode("Element for face identification is not a face");
1122     return NULL;
1123   }
1124
1125   //Compute the Face value
1126   try {
1127 #if OCC_VERSION_LARGE > 0x06010000
1128     OCC_CATCH_SIGNALS;
1129 #endif
1130     TopoDS_Shape aShape;
1131
1132     GEOMImpl_Block6Explorer aBlockTool;
1133     aBlockTool.InitByBlockAndFace(aBlockOrComp, anArg);
1134     aShape = aBlockTool.GetFace(2);
1135
1136     TopTools_IndexedMapOfShape anIndices;
1137     TopExp::MapShapes(aBlockOrComp, anIndices);
1138     Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
1139     anArray->SetValue(1, anIndices.FindIndex(aShape));
1140     aResult = GetEngine()->AddSubShape(theShape, anArray);
1141   }
1142   catch (Standard_Failure) {
1143     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1144     SetErrorCode(aFail->GetMessageString());
1145     return NULL;
1146   }
1147
1148   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1149
1150   //Make a Python command
1151   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetOppositeFace("
1152                                << theShape << ", " << theFace << ")";
1153
1154   SetErrorCode(OK);
1155   return aResult;
1156 }
1157
1158 //=============================================================================
1159 /*!
1160  *  GetFaceNearPoint
1161  */
1162 //=============================================================================
1163 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceNearPoint
1164                                                 (Handle(GEOM_Object) theShape,
1165                                                  Handle(GEOM_Object) thePoint)
1166 {
1167   SetErrorCode(KO);
1168
1169   //New object
1170   Handle(GEOM_Object) aResult;
1171
1172   // Arguments
1173   if (theShape.IsNull() || thePoint.IsNull()) return NULL;
1174
1175   TopoDS_Shape aBlockOrComp = theShape->GetValue();
1176   if (aBlockOrComp.IsNull()) {
1177     SetErrorCode("Block or compound is null");
1178     return NULL;
1179   }
1180
1181   TopoDS_Shape anArg = thePoint->GetValue();
1182   if (anArg.IsNull()) {
1183     SetErrorCode("Null shape is given as argument");
1184     return NULL;
1185   }
1186   if (anArg.ShapeType() != TopAbs_VERTEX) {
1187     SetErrorCode("Element for face identification is not a vertex");
1188     return NULL;
1189   }
1190
1191   //Compute the Face value
1192   try {
1193 #if OCC_VERSION_LARGE > 0x06010000
1194     OCC_CATCH_SIGNALS;
1195 #endif
1196     TopoDS_Shape aShape;
1197
1198     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
1199     gp_Pnt aPnt = BRep_Tool::Pnt(aVert);
1200     Standard_Real PX, PY, PZ;
1201     aPnt.Coord(PX, PY, PZ);
1202
1203     // 1. Classify the point relatively each face
1204     Standard_Integer nearest = 2, nbFound = 0;
1205     TopTools_DataMapOfShapeInteger mapShapeDist;
1206     TopExp_Explorer exp (aBlockOrComp, TopAbs_FACE);
1207     for (; exp.More(); exp.Next()) {
1208       TopoDS_Shape aFace = exp.Current();
1209
1210       if (!mapShapeDist.IsBound(aFace)) {
1211         Standard_Integer aDistance = 2;
1212
1213         // 1.a. Classify relatively Surface
1214         Handle(Geom_Surface) aSurf = BRep_Tool::Surface(TopoDS::Face(aFace));
1215         Handle(ShapeAnalysis_Surface) aSurfAna = new ShapeAnalysis_Surface (aSurf);
1216         gp_Pnt2d p2dOnSurf = aSurfAna->ValueOfUV(aPnt, Precision::Confusion());
1217         gp_Pnt p3dOnSurf = aSurfAna->Value(p2dOnSurf);
1218         Standard_Real aDist = p3dOnSurf.Distance(aPnt);
1219         if (aDist > Precision::Confusion()) {
1220           // OUT of Surface
1221           aDistance = 1;
1222         } else {
1223           // 1.b. Classify relatively the face itself
1224           BRepClass_FaceClassifier FC (TopoDS::Face(aFace), p2dOnSurf, Precision::Confusion());
1225           if (FC.State() == TopAbs_IN) {
1226             aDistance = -1;
1227           } else if (FC.State() == TopAbs_ON) {
1228             aDistance = 0;
1229           } else { // OUT
1230             aDistance = 1;
1231           }
1232         }
1233
1234         if (aDistance < nearest) {
1235           nearest = aDistance;
1236           aShape = aFace;
1237           nbFound = 1;
1238
1239           // A first found face, containing the point inside, will be returned.
1240           // It is the solution, if there are no
1241           // coincident or intersecting faces in the compound.
1242           if (nearest == -1) break;
1243
1244         } else if (aDistance == nearest) {
1245           nbFound++;
1246         } else {
1247         }
1248
1249         mapShapeDist.Bind(aFace, aDistance);
1250       } // if (!mapShapeDist.IsBound(aFace))
1251     }
1252
1253     // 2. Define face, containing the point or having minimum distance to it
1254     if (nbFound > 1) {
1255       if (nearest == 0) {
1256         // The point is on boundary of some faces and there are
1257         // no faces, having the point inside
1258         SetErrorCode("Multiple faces near the given point are found");
1259         return NULL;
1260
1261       } else if (nearest == 1) {
1262         // The point is outside some faces and there are
1263         // no faces, having the point inside or on boundary.
1264         // We will get a nearest face
1265         Standard_Real bigReal = RealLast();
1266         Standard_Real minDist = bigReal;
1267         TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
1268         for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
1269           if (mapShapeDistIter.Value() == 1) {
1270             TopoDS_Shape aFace = mapShapeDistIter.Key();
1271             Standard_Real aDist = bigReal;
1272
1273             // 2.a. Fast check of distance - if point projection on surface is on face
1274             Handle(Geom_Surface) aSurf = BRep_Tool::Surface(TopoDS::Face(aFace));
1275             Handle(ShapeAnalysis_Surface) aSurfAna = new ShapeAnalysis_Surface (aSurf);
1276             gp_Pnt2d p2dOnSurf = aSurfAna->ValueOfUV(aPnt, Precision::Confusion());
1277             gp_Pnt p3dOnSurf = aSurfAna->Value(p2dOnSurf);
1278             aDist = p3dOnSurf.Distance(aPnt);
1279
1280             BRepClass_FaceClassifier FC (TopoDS::Face(aFace), p2dOnSurf, Precision::Confusion());
1281             if (FC.State() == TopAbs_OUT) {
1282               if (aDist < minDist) {
1283                 // 2.b. Slow check - if point projection on surface is outside of face
1284                 BRepExtrema_DistShapeShape aDistTool (aVert, aFace);
1285                 if (!aDistTool.IsDone()) {
1286                   SetErrorCode("Can not find a distance from the given point to one of faces");
1287                   return NULL;
1288                 }
1289                 aDist = aDistTool.Value();
1290               } else {
1291                 aDist = bigReal;
1292               }
1293             }
1294
1295             if (aDist < minDist) {
1296               minDist = aDist;
1297               aShape = aFace;
1298             }
1299           }
1300         }
1301       } else { // nearest == -1
1302 //        // The point is inside some faces.
1303 //        // We will get a face with nearest center
1304 //        Standard_Real minDist = RealLast();
1305 //        TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
1306 //        for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
1307 //          if (mapShapeDistIter.Value() == -1) {
1308 //            TopoDS_Shape aFace = mapShapeDistIter.Key();
1309 //            GProp_GProps aSystem;
1310 //            BRepGProp::SurfaceProperties(aFace, aSystem);
1311 //            gp_Pnt aCenterMass = aSystem.CentreOfMass();
1312 //
1313 //            Standard_Real aDist = aCenterMass.Distance(aPnt);
1314 //            if (aDist < minDist) {
1315 //              minDist = aDist;
1316 //              aShape = aFace;
1317 //            }
1318 //          }
1319 //        }
1320       }
1321     } // if (nbFound > 1)
1322
1323     if (nbFound == 0) {
1324       SetErrorCode("There are no faces near the given point");
1325       return NULL;
1326     } else {
1327       TopTools_IndexedMapOfShape anIndices;
1328       TopExp::MapShapes(aBlockOrComp, anIndices);
1329       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
1330       anArray->SetValue(1, anIndices.FindIndex(aShape));
1331       aResult = GetEngine()->AddSubShape(theShape, anArray);
1332     }
1333   }
1334   catch (Standard_Failure) {
1335     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1336     SetErrorCode(aFail->GetMessageString());
1337     return NULL;
1338   }
1339
1340   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1341
1342   //Make a Python command
1343   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceNearPoint("
1344                                << theShape << ", " << thePoint << ")";
1345
1346   SetErrorCode(OK);
1347   return aResult;
1348 }
1349
1350 //=============================================================================
1351 /*!
1352  *  GetFaceByNormale
1353  */
1354 //=============================================================================
1355 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByNormale
1356                                                 (Handle(GEOM_Object) theShape,
1357                                                  Handle(GEOM_Object) theVector)
1358 {
1359   SetErrorCode(KO);
1360
1361   //New object
1362   Handle(GEOM_Object) aResult;
1363
1364   // Arguments
1365   if (theShape.IsNull() || theVector.IsNull()) return NULL;
1366
1367   TopoDS_Shape aBlockOrComp = theShape->GetValue();
1368   if (aBlockOrComp.IsNull()) {
1369     SetErrorCode("Block or compound is null");
1370     return NULL;
1371   }
1372
1373   TopoDS_Shape anArg = theVector->GetValue();
1374   if (anArg.IsNull()) {
1375     SetErrorCode("Null shape is given as argument");
1376     return NULL;
1377   }
1378   if (anArg.ShapeType() != TopAbs_EDGE) {
1379     SetErrorCode("Element for normale identification is not an edge");
1380     return NULL;
1381   }
1382
1383   //Compute the Face value
1384   try {
1385 #if OCC_VERSION_LARGE > 0x06010000
1386     OCC_CATCH_SIGNALS;
1387 #endif
1388     TopoDS_Shape aShape;
1389
1390     TopoDS_Edge anEdge = TopoDS::Edge(anArg);
1391     TopoDS_Vertex V1, V2;
1392     TopExp::Vertices(anEdge, V1, V2, Standard_True);
1393     gp_Pnt P1 = BRep_Tool::Pnt(V1);
1394     gp_Pnt P2 = BRep_Tool::Pnt(V2);
1395     gp_Vec aVec (P1, P2);
1396     if (aVec.Magnitude() < Precision::Confusion()) {
1397       SetErrorCode("Vector with null magnitude is given");
1398       return NULL;
1399     }
1400
1401     Standard_Real minAngle = RealLast();
1402     TopTools_MapOfShape mapShape;
1403     TopExp_Explorer exp (aBlockOrComp, TopAbs_FACE);
1404     for (; exp.More(); exp.Next()) {
1405       if (mapShape.Add(exp.Current())) {
1406         TopoDS_Face aFace = TopoDS::Face(exp.Current());
1407         BRepAdaptor_Surface SF (aFace);
1408
1409         Standard_Real u, v, x;
1410
1411         // find a point on the surface to get normal direction in
1412         u = SF.FirstUParameter();
1413         x = SF.LastUParameter();
1414         if (Precision::IsInfinite(u)) {
1415           u =  (Precision::IsInfinite(x)) ? 0. : x;
1416         } else if (!Precision::IsInfinite(x)) {
1417           u = (u+x) / 2.;
1418         }
1419
1420         v = SF.FirstVParameter();
1421         x = SF.LastVParameter();
1422         if (Precision::IsInfinite(v)) {
1423           v =  (Precision::IsInfinite(x)) ? 0. : x;
1424         } else if (!Precision::IsInfinite(x)) {
1425           v = (v+x) / 2.;
1426         }
1427
1428         // compute the normal direction
1429         gp_Vec Vec1,Vec2;
1430         SF.D1(u,v,P1,Vec1,Vec2);
1431         gp_Vec V = Vec1.Crossed(Vec2);
1432         x = V.Magnitude();
1433         if (V.Magnitude() < Precision::Confusion()) {
1434           SetErrorCode("Normal vector of a face has null magnitude");
1435           return NULL;
1436         }
1437
1438         // consider the face orientation
1439         if (aFace.Orientation() == TopAbs_REVERSED ||
1440             aFace.Orientation() == TopAbs_INTERNAL) {
1441           V = - V;
1442         }
1443
1444         // compute the angle and compare with the minimal one
1445         Standard_Real anAngle = aVec.Angle(V);
1446         if (anAngle < minAngle) {
1447           minAngle = anAngle;
1448           aShape = aFace;
1449         }
1450       }
1451     }
1452
1453     if (aShape.IsNull()) {
1454       SetErrorCode("Failed to find a face by the given normale");
1455       return NULL;
1456     } else {
1457       TopTools_IndexedMapOfShape anIndices;
1458       TopExp::MapShapes(aBlockOrComp, anIndices);
1459       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
1460       anArray->SetValue(1, anIndices.FindIndex(aShape));
1461       aResult = GetEngine()->AddSubShape(theShape, anArray);
1462     }
1463   }
1464   catch (Standard_Failure) {
1465     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1466     SetErrorCode(aFail->GetMessageString());
1467     return NULL;
1468   }
1469
1470   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1471
1472   //Make a Python command
1473   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceByNormale("
1474     << theShape << ", " << theVector << ")";
1475
1476   SetErrorCode(OK);
1477   return aResult;
1478 }
1479
1480 //=============================================================================
1481 /*!
1482  *  GetShapesNearPoint
1483  */
1484 //=============================================================================
1485 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetShapesNearPoint
1486                                          (Handle(GEOM_Object)    theShape,
1487                                           Handle(GEOM_Object)    thePoint,
1488                                           const Standard_Integer theShapeType,
1489                                           const Standard_Real    theConstTolerance)
1490 {
1491   SetErrorCode(KO);
1492
1493   // New object
1494   Handle(GEOM_Object) aResult;
1495
1496   // Arguments
1497   if (theShape.IsNull() || thePoint.IsNull()) return NULL;
1498
1499   TopoDS_Shape aBlockOrComp = theShape->GetValue();
1500   if (aBlockOrComp.IsNull()) {
1501     SetErrorCode("Block or compound is null");
1502     return NULL;
1503   }
1504
1505   TopoDS_Shape anArg = thePoint->GetValue();
1506   if (anArg.IsNull()) {
1507     SetErrorCode("Null shape is given as argument");
1508     return NULL;
1509   }
1510   if (anArg.ShapeType() != TopAbs_VERTEX) {
1511     SetErrorCode("Element for face identification is not a vertex");
1512     return NULL;
1513   }
1514
1515   if (theShapeType < TopAbs_SOLID || TopAbs_VERTEX < theShapeType) {
1516     SetErrorCode("Invalid type of result is requested");
1517     return NULL;
1518   }
1519
1520   Standard_Real theTolerance = theConstTolerance;
1521   if (theTolerance < Precision::Confusion()) {
1522     theTolerance = Precision::Confusion();
1523   }
1524
1525   // Compute the result
1526   try {
1527 #if OCC_VERSION_LARGE > 0x06010000
1528     OCC_CATCH_SIGNALS;
1529 #endif
1530     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
1531
1532     TopTools_MapOfShape mapShape;
1533     Standard_Integer nbEdges = 0;
1534     TopExp_Explorer exp (aBlockOrComp, TopAbs_ShapeEnum(theShapeType));
1535     for (; exp.More(); exp.Next()) {
1536       if (mapShape.Add(exp.Current())) {
1537         nbEdges++;
1538       }
1539     }
1540
1541     if (nbEdges == 0) {
1542       SetErrorCode("Given shape contains no sub-shapes of requested type");
1543       return NULL;
1544     }
1545
1546     // Calculate distances and find min
1547     mapShape.Clear();
1548     Standard_Integer ind = 1;
1549     Standard_Real aMinDist = RealLast();
1550     TopTools_Array1OfShape anEdges (1, nbEdges);
1551     TColStd_Array1OfReal aDistances (1, nbEdges);
1552     for (exp.Init(aBlockOrComp, TopAbs_ShapeEnum(theShapeType)); exp.More(); exp.Next()) {
1553       if (mapShape.Add(exp.Current())) {
1554         TopoDS_Shape anEdge = exp.Current();
1555         anEdges(ind) = anEdge;
1556
1557         BRepExtrema_DistShapeShape aDistTool (aVert, anEdges(ind));
1558         if (!aDistTool.IsDone()) {
1559           SetErrorCode("Can not find a distance from the given point to one of sub-shapes");
1560           return NULL;
1561         }
1562         aDistances(ind) = aDistTool.Value();
1563         if (aDistances(ind) < aMinDist) {
1564           aMinDist = aDistances(ind);
1565         }
1566         ind++;
1567       }
1568     }
1569
1570     if (aMinDist < RealLast()) {
1571       // Collect sub-shapes with distance < (aMinDist + theTolerance)
1572       int nbSubShapes = 0;
1573       TopTools_Array1OfShape aNearShapes (1, nbEdges);
1574       for (ind = 1; ind <= nbEdges; ind++) {
1575         if (aDistances(ind) < aMinDist + theTolerance) {
1576           nbSubShapes++;
1577           aNearShapes(nbSubShapes) = anEdges(ind);
1578         }
1579       }
1580
1581       // Add sub-shape
1582       TopTools_IndexedMapOfShape anIndices;
1583       TopExp::MapShapes(aBlockOrComp, anIndices);
1584       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger (1, nbSubShapes);
1585       for (ind = 1; ind <= nbSubShapes; ind++) {
1586         anArray->SetValue(ind, anIndices.FindIndex(aNearShapes(ind)));
1587       }
1588       aResult = GetEngine()->AddSubShape(theShape, anArray);
1589     }
1590   }
1591   catch (Standard_Failure) {
1592     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1593     SetErrorCode(aFail->GetMessageString());
1594     return NULL;
1595   }
1596
1597   if (aResult.IsNull())
1598     return NULL;
1599
1600   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1601
1602   //Make a Python command
1603   GEOM::TPythonDump(aFunction)
1604     << aResult << " = geompy.GetShapesNearPoint(" << theShape << ", " << thePoint
1605     << ", " << TopAbs_ShapeEnum(theShapeType) << ", " << theTolerance << ")";
1606
1607   SetErrorCode(OK);
1608   return aResult;
1609 }
1610
1611 //=============================================================================
1612 /*!
1613  *  IsCompoundOfBlocks
1614  */
1615 //=============================================================================
1616 Standard_Boolean GEOMImpl_IBlocksOperations::IsCompoundOfBlocks
1617                                                 (Handle(GEOM_Object)    theCompound,
1618                                                  const Standard_Integer theMinNbFaces,
1619                                                  const Standard_Integer theMaxNbFaces,
1620                                                  Standard_Integer&      theNbBlocks)
1621 {
1622   SetErrorCode(KO);
1623   Standard_Boolean isCompOfBlocks = Standard_False;
1624   theNbBlocks = 0;
1625
1626   if (theCompound.IsNull()) return isCompOfBlocks;
1627   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
1628
1629   //Check
1630   isCompOfBlocks = Standard_True;
1631   try {
1632 #if OCC_VERSION_LARGE > 0x06010000
1633     OCC_CATCH_SIGNALS;
1634 #endif
1635     TopTools_MapOfShape mapShape;
1636     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
1637     for (; exp.More(); exp.Next()) {
1638       if (mapShape.Add(exp.Current())) {
1639         TopoDS_Shape aSolid = exp.Current();
1640
1641         TopTools_MapOfShape mapFaces;
1642         TopExp_Explorer expF (aSolid, TopAbs_FACE);
1643         Standard_Integer nbFaces = 0;
1644         for (; expF.More(); expF.Next()) {
1645           if (mapFaces.Add(expF.Current())) {
1646             nbFaces++;
1647             if (nbFaces > theMaxNbFaces) {
1648               isCompOfBlocks = Standard_False;
1649               break;
1650             }
1651           }
1652         }
1653         if (nbFaces < theMinNbFaces || theMaxNbFaces < nbFaces) {
1654           isCompOfBlocks = Standard_False;
1655         } else {
1656           theNbBlocks++;
1657         }
1658       }
1659     }
1660   }
1661   catch (Standard_Failure) {
1662     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1663     SetErrorCode(aFail->GetMessageString());
1664     return isCompOfBlocks;
1665   }
1666
1667   SetErrorCode(OK);
1668   return isCompOfBlocks;
1669 }
1670
1671 //=============================================================================
1672 /*!
1673  *  Set of functions, used by CheckCompoundOfBlocks() method
1674  */
1675 //=============================================================================
1676 void GEOMImpl_IBlocksOperations::AddBlocksFrom (const TopoDS_Shape&   theShape,
1677                                                 TopTools_ListOfShape& BLO,
1678                                                 TopTools_ListOfShape& NOT,
1679                                                 TopTools_ListOfShape& EXT,
1680                                                 TopTools_ListOfShape& NOQ)
1681 {
1682   TopAbs_ShapeEnum aType = theShape.ShapeType();
1683   switch (aType) {
1684   case TopAbs_COMPOUND:
1685   case TopAbs_COMPSOLID:
1686     {
1687       TopoDS_Iterator It (theShape);
1688       for (; It.More(); It.Next()) {
1689         AddBlocksFrom(It.Value(), BLO, NOT, EXT, NOQ);
1690       }
1691     }
1692     break;
1693   case TopAbs_SOLID:
1694     {
1695       // Check, if there are seam or degenerated edges
1696       BlockFix_CheckTool aTool;
1697       aTool.SetShape(theShape);
1698       aTool.Perform();
1699       if (aTool.NbPossibleBlocks() > 0) {
1700         EXT.Append(theShape);
1701       } else {
1702         // Count faces and edges in each face to recognize blocks
1703         TopTools_MapOfShape mapFaces;
1704         Standard_Integer nbFaces = 0;
1705         Standard_Boolean hasNonQuadr = Standard_False;
1706         TopExp_Explorer expF (theShape, TopAbs_FACE);
1707
1708         for (; expF.More(); expF.Next()) {
1709           if (mapFaces.Add(expF.Current())) {
1710             nbFaces++;
1711             //0021483//if (nbFaces > 6) break;
1712
1713             // get wire
1714             TopoDS_Shape aF = expF.Current();
1715             TopExp_Explorer wires (aF, TopAbs_WIRE);
1716             if (!wires.More()) {
1717               // no wire in the face
1718               hasNonQuadr = Standard_True;
1719               NOQ.Append(aF);//0021483
1720               //0021483//break;
1721               continue;
1722             }
1723             TopoDS_Shape aWire = wires.Current();
1724             wires.Next();
1725             if (wires.More()) {
1726               // multiple wires in the face
1727               hasNonQuadr = Standard_True;
1728               NOQ.Append(aF);//0021483
1729               //0021483//break;
1730               continue;
1731             }
1732
1733             // Check number of edges in the face
1734             Standard_Integer nbEdges = 0;
1735             TopTools_MapOfShape mapEdges;
1736             TopExp_Explorer expW (aWire, TopAbs_EDGE);
1737             for (; expW.More(); expW.Next()) {
1738               if (mapEdges.Add(expW.Current())) {
1739                 nbEdges++;
1740                 if (nbEdges > 4) break;
1741               }
1742             }
1743             if (nbEdges != 4) {
1744               hasNonQuadr = Standard_True;
1745               NOQ.Append(aF);//0021483
1746             }
1747           }
1748         }
1749
1750         if (nbFaces == 6 && !hasNonQuadr) {
1751           BLO.Append(theShape);
1752         } else {
1753           NOT.Append(theShape);
1754         }
1755       }
1756     }
1757     break;
1758   case TopAbs_SHELL: //0021483
1759   case TopAbs_FACE: //0021483
1760     {
1761       // Count edges in each face
1762       TopTools_MapOfShape mapFaces;
1763       TopExp_Explorer expF (theShape, TopAbs_FACE);
1764       for (; expF.More(); expF.Next()) {
1765         if (mapFaces.Add(expF.Current())) {
1766           // get wire
1767           TopoDS_Shape aF = expF.Current();
1768           TopExp_Explorer wires (aF, TopAbs_WIRE);
1769           if (!wires.More()) {
1770             // no wire in the face
1771             NOQ.Append(aF);//0021483
1772             continue;
1773           }
1774           TopoDS_Shape aWire = wires.Current();
1775           wires.Next();
1776           if (wires.More()) {
1777             // multiple wires in the face
1778             NOQ.Append(aF);//0021483
1779             continue;
1780           }
1781
1782           // Check number of edges in the face
1783           Standard_Integer nbEdges = 0;
1784           TopTools_MapOfShape mapEdges;
1785           TopExp_Explorer expW (aWire, TopAbs_EDGE);
1786           for (; expW.More(); expW.Next()) {
1787             if (mapEdges.Add(expW.Current())) {
1788               nbEdges++;
1789               if (nbEdges > 4) break;
1790             }
1791           }
1792           if (nbEdges != 4) {
1793             NOQ.Append(aF);//0021483
1794           }
1795         }
1796       }
1797     }
1798     break;
1799   default:
1800     NOT.Append(theShape);
1801   }
1802 }
1803
1804 void AddBlocksFromOld (const TopoDS_Shape&   theShape,
1805                        TopTools_ListOfShape& BLO,
1806                        TopTools_ListOfShape& NOT,
1807                        TopTools_ListOfShape& DEG,
1808                        TopTools_ListOfShape& SEA)
1809 {
1810   TopAbs_ShapeEnum aType = theShape.ShapeType();
1811   switch (aType) {
1812   case TopAbs_COMPOUND:
1813   case TopAbs_COMPSOLID:
1814     {
1815       TopoDS_Iterator It (theShape);
1816       for (; It.More(); It.Next()) {
1817         AddBlocksFromOld(It.Value(), BLO, NOT, DEG, SEA);
1818       }
1819     }
1820     break;
1821   case TopAbs_SOLID:
1822     {
1823       TopTools_MapOfShape mapFaces;
1824       TopExp_Explorer expF (theShape, TopAbs_FACE);
1825       Standard_Integer nbFaces = 0;
1826       Standard_Boolean hasNonQuadr = Standard_False;
1827       Standard_Boolean hasDegenerated = Standard_False;
1828       Standard_Boolean hasSeam = Standard_False;
1829       for (; expF.More(); expF.Next()) {
1830         if (mapFaces.Add(expF.Current())) {
1831           nbFaces++;
1832           if (nbFaces > 6) break;
1833
1834           // Check number of edges in the face
1835           Standard_Integer nbEdges = 0;
1836           TopTools_MapOfShape mapEdges;
1837
1838           // get wire
1839           TopoDS_Shape aF = expF.Current();
1840           TopExp_Explorer wires (aF, TopAbs_WIRE);
1841           if (!wires.More()) {
1842             // no wire in the face
1843             hasNonQuadr = Standard_True;
1844             break;
1845           }
1846           TopoDS_Shape aWire = wires.Current();
1847           wires.Next();
1848           if (wires.More()) {
1849             // multiple wires in the face
1850             hasNonQuadr = Standard_True;
1851             break;
1852           }
1853
1854           // iterate on wire
1855           BRepTools_WireExplorer aWE (TopoDS::Wire(aWire), TopoDS::Face(aF));
1856           for (; aWE.More(); aWE.Next(), nbEdges++) {
1857             if (BRep_Tool::Degenerated(aWE.Current())) {
1858               // degenerated edge found
1859               hasDegenerated = Standard_True;
1860 //              break;
1861             }
1862             if (mapEdges.Contains(aWE.Current())) {
1863               // seam edge found
1864               hasSeam = Standard_True;
1865 //              break;
1866             }
1867             mapEdges.Add(aWE.Current());
1868           }
1869           if (nbEdges != 4) {
1870             hasNonQuadr = Standard_True;
1871           }
1872         }
1873       }
1874       if (nbFaces == 6) {
1875         if (hasDegenerated || hasSeam) {
1876           if (hasDegenerated) {
1877             DEG.Append(theShape);
1878           }
1879           if (hasSeam) {
1880             SEA.Append(theShape);
1881           }
1882         } else if (hasNonQuadr) {
1883           NOT.Append(theShape);
1884         } else {
1885           BLO.Append(theShape);
1886         }
1887       } else {
1888         NOT.Append(theShape);
1889       }
1890     }
1891     break;
1892   default:
1893     NOT.Append(theShape);
1894   }
1895 }
1896
1897 #define REL_NOT_CONNECTED 0
1898 #define REL_OK            1
1899 #define REL_NOT_GLUED     2
1900 #define REL_COLLISION_VV  3
1901 #define REL_COLLISION_FF  4
1902 #define REL_COLLISION_EE  5
1903 #define REL_UNKNOWN       6
1904
1905 Standard_Integer BlocksRelation (const TopoDS_Shape& theBlock1,
1906                                  const TopoDS_Shape& theBlock2)
1907 {
1908   // Compare bounding boxes before calling BRepExtrema_DistShapeShape
1909   Standard_Real Xmin1, Ymin1, Zmin1, Xmax1, Ymax1, Zmax1;
1910   Standard_Real Xmin2, Ymin2, Zmin2, Xmax2, Ymax2, Zmax2;
1911   Bnd_Box B1, B2;
1912   BRepBndLib::Add(theBlock1, B1);
1913   BRepBndLib::Add(theBlock2, B2);
1914   B1.Get(Xmin1, Ymin1, Zmin1, Xmax1, Ymax1, Zmax1);
1915   B2.Get(Xmin2, Ymin2, Zmin2, Xmax2, Ymax2, Zmax2);
1916   if (Xmax2 < Xmin1 || Xmax1 < Xmin2 ||
1917       Ymax2 < Ymin1 || Ymax1 < Ymin2 ||
1918       Zmax2 < Zmin1 || Zmax1 < Zmin2) {
1919     return REL_NOT_CONNECTED;
1920   }
1921
1922   BRepExtrema_DistShapeShape dst (theBlock1, theBlock2);
1923   if (!dst.IsDone()) {
1924     return REL_UNKNOWN;
1925   }
1926
1927   if (dst.Value() > Precision::Confusion()) {
1928     return REL_NOT_CONNECTED;
1929   }
1930
1931   if (dst.InnerSolution()) {
1932     return REL_COLLISION_VV;
1933   }
1934
1935   Standard_Integer nbSol = dst.NbSolution();
1936   Standard_Integer relation = REL_OK;
1937   Standard_Integer nbVerts = 0;
1938   Standard_Integer nbEdges = 0;
1939   Standard_Integer sol = 1;
1940   for (; sol <= nbSol; sol++) {
1941     BRepExtrema_SupportType supp1 = dst.SupportTypeShape1(sol);
1942     BRepExtrema_SupportType supp2 = dst.SupportTypeShape2(sol);
1943     if (supp1 == BRepExtrema_IsVertex && supp2 == BRepExtrema_IsVertex) {
1944       nbVerts++;
1945     } else if (supp1 == BRepExtrema_IsInFace || supp2 == BRepExtrema_IsInFace) {
1946       return REL_COLLISION_FF;
1947     } else if (supp1 == BRepExtrema_IsOnEdge && supp2 == BRepExtrema_IsOnEdge) {
1948       nbEdges++;
1949     } else if ((supp1 == BRepExtrema_IsOnEdge && supp2 == BRepExtrema_IsVertex) ||
1950                (supp2 == BRepExtrema_IsOnEdge && supp1 == BRepExtrema_IsVertex)) {
1951       relation = REL_COLLISION_EE;
1952     } else {
1953     }
1954   }
1955
1956   if (relation != REL_OK) {
1957     return relation;
1958   }
1959
1960   TColStd_Array1OfInteger vertSol (1, nbVerts);
1961   TopTools_Array1OfShape V1 (1, nbVerts);
1962   TopTools_Array1OfShape V2 (1, nbVerts);
1963   Standard_Integer ivs = 0;
1964   for (sol = 1; sol <= nbSol; sol++) {
1965     if (dst.SupportTypeShape1(sol) == BRepExtrema_IsVertex &&
1966         dst.SupportTypeShape2(sol) == BRepExtrema_IsVertex) {
1967       TopoDS_Vertex Vcur = TopoDS::Vertex(dst.SupportOnShape1(sol));
1968       // Check, that this vertex is far enough from other solution vertices.
1969       Standard_Integer ii = 1;
1970       for (; ii <= ivs; ii++) {
1971         if (BRepTools::Compare(TopoDS::Vertex(V1(ii)), Vcur)) {
1972           continue;
1973         }
1974       }
1975       ivs++;
1976       vertSol(ivs) = sol;
1977       V1(ivs) = Vcur;
1978       V2(ivs) = dst.SupportOnShape2(sol);
1979     }
1980   }
1981
1982   // As we deal only with quadrangles,
1983   // 2, 3 or 4 vertex solutions can be found.
1984   if (ivs <= 1) {
1985     if (nbEdges > 0) {
1986       return REL_COLLISION_FF;
1987     }
1988     return REL_NOT_CONNECTED;
1989   }
1990   if (ivs > 4) {
1991     return REL_UNKNOWN;
1992   }
1993
1994   // Check sharing of coincident entities.
1995   if (ivs == 2 || ivs == 3) {
1996     // Map vertices and edges of the blocks
1997     TopTools_IndexedDataMapOfShapeListOfShape MVE1, MVE2;
1998     GEOMImpl_Block6Explorer::MapShapesAndAncestors
1999       (theBlock1, TopAbs_VERTEX, TopAbs_EDGE, MVE1);
2000     GEOMImpl_Block6Explorer::MapShapesAndAncestors
2001       (theBlock2, TopAbs_VERTEX, TopAbs_EDGE, MVE2);
2002
2003     if (ivs == 2) {
2004       // Find common edge
2005       TopoDS_Shape anEdge1, anEdge2;
2006       GEOMImpl_Block6Explorer::FindEdge(anEdge1, V1(1), V1(2), MVE1);
2007       if (anEdge1.IsNull()) return REL_UNKNOWN;
2008
2009       GEOMImpl_Block6Explorer::FindEdge(anEdge2, V2(1), V2(2), MVE2);
2010       if (anEdge2.IsNull()) return REL_UNKNOWN;
2011
2012       if (!anEdge1.IsSame(anEdge2)) return REL_NOT_GLUED;
2013
2014     } else { // ivs == 3
2015       // Find common edges
2016       Standard_Integer e1_v1 = 1;
2017       Standard_Integer e1_v2 = 2;
2018       Standard_Integer e2_v1 = 3;
2019       Standard_Integer e2_v2 = 1;
2020
2021       TopoDS_Shape anEdge11, anEdge12;
2022       GEOMImpl_Block6Explorer::FindEdge(anEdge11, V1(e1_v1), V1(e1_v2), MVE1);
2023       if (anEdge11.IsNull()) {
2024         e1_v2 = 3;
2025         e2_v1 = 2;
2026         GEOMImpl_Block6Explorer::FindEdge(anEdge11, V1(e1_v1), V1(e1_v2), MVE1);
2027         if (anEdge11.IsNull()) return REL_UNKNOWN;
2028       }
2029       GEOMImpl_Block6Explorer::FindEdge(anEdge12, V1(e2_v1), V1(e2_v2), MVE1);
2030       if (anEdge12.IsNull()) {
2031         e2_v2 = 5 - e2_v1;
2032         GEOMImpl_Block6Explorer::FindEdge(anEdge12, V1(e2_v1), V1(e2_v2), MVE1);
2033         if (anEdge12.IsNull()) return REL_UNKNOWN;
2034       }
2035
2036       TopoDS_Shape anEdge21, anEdge22;
2037       GEOMImpl_Block6Explorer::FindEdge(anEdge21, V2(e1_v1), V2(e1_v2), MVE2);
2038       if (anEdge21.IsNull()) return REL_UNKNOWN;
2039       GEOMImpl_Block6Explorer::FindEdge(anEdge22, V2(e2_v1), V2(e2_v2), MVE2);
2040       if (anEdge22.IsNull()) return REL_UNKNOWN;
2041
2042       // Check of edges coincidence (with some precision) have to be done here
2043       // if (!anEdge11.IsEqual(anEdge21)) return REL_UNKNOWN;
2044       // if (!anEdge12.IsEqual(anEdge22)) return REL_UNKNOWN;
2045
2046       // Check of edges sharing
2047       if (!anEdge11.IsSame(anEdge21)) return REL_NOT_GLUED;
2048       if (!anEdge12.IsSame(anEdge22)) return REL_NOT_GLUED;
2049     }
2050   }
2051
2052   if (ivs == 4) {
2053     // Map vertices and faces of the blocks
2054     TopTools_IndexedDataMapOfShapeListOfShape MVF1, MVF2;
2055     GEOMImpl_Block6Explorer::MapShapesAndAncestors
2056       (theBlock1, TopAbs_VERTEX, TopAbs_FACE, MVF1);
2057     GEOMImpl_Block6Explorer::MapShapesAndAncestors
2058       (theBlock2, TopAbs_VERTEX, TopAbs_FACE, MVF2);
2059
2060     TopoDS_Shape aFace1, aFace2;
2061     GEOMImpl_Block6Explorer::FindFace(aFace1, V1(1), V1(2), V1(3), V1(4), MVF1);
2062     if (aFace1.IsNull()) return REL_UNKNOWN;
2063     GEOMImpl_Block6Explorer::FindFace(aFace2, V2(1), V2(2), V2(3), V2(4), MVF2);
2064     if (aFace2.IsNull()) return REL_UNKNOWN;
2065
2066     // Check of faces coincidence (with some precision) have to be done here
2067     // if (!aFace1.IsEqual(aFace2)) return REL_UNKNOWN;
2068
2069     // Check of faces sharing
2070     if (!aFace1.IsSame(aFace2)) return REL_NOT_GLUED;
2071   }
2072
2073   return REL_OK;
2074 }
2075
2076 void FindConnected (const Standard_Integer         theBlockIndex,
2077                     const TColStd_Array2OfInteger& theRelations,
2078                     TColStd_MapOfInteger&          theProcessedMap,
2079                     TColStd_MapOfInteger&          theConnectedMap)
2080 {
2081   theConnectedMap.Add(theBlockIndex);
2082   theProcessedMap.Add(theBlockIndex);
2083
2084   Standard_Integer nbBlocks = theRelations.ColLength();
2085   Standard_Integer col = 1;
2086   for (; col <= nbBlocks; col++) {
2087     if (theRelations(theBlockIndex, col) == REL_OK ||
2088         theRelations(theBlockIndex, col) == REL_NOT_GLUED) {
2089       if (!theProcessedMap.Contains(col)) {
2090         FindConnected(col, theRelations, theProcessedMap, theConnectedMap);
2091       }
2092     }
2093   }
2094 }
2095
2096 Standard_Boolean HasAnyConnection (const Standard_Integer         theBlockIndex,
2097                                    const TColStd_MapOfInteger&    theWith,
2098                                    const TColStd_Array2OfInteger& theRelations,
2099                                    TColStd_MapOfInteger&          theProcessedMap)
2100 {
2101   theProcessedMap.Add(theBlockIndex);
2102
2103   Standard_Integer nbBlocks = theRelations.ColLength();
2104   Standard_Integer col = 1;
2105   for (; col <= nbBlocks; col++) {
2106     if (theRelations(theBlockIndex, col) != REL_NOT_CONNECTED) {
2107       if (!theProcessedMap.Contains(col)) {
2108         if (theWith.Contains(col))
2109           return Standard_True;
2110         if (HasAnyConnection(col, theWith, theRelations, theProcessedMap))
2111           return Standard_True;
2112       }
2113     }
2114   }
2115
2116   return Standard_False;
2117 }
2118
2119 //=============================================================================
2120 /*!
2121  *  CheckCompoundOfBlocksOld
2122  */
2123 //=============================================================================
2124 Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocksOld
2125                                                 (Handle(GEOM_Object) theCompound,
2126                                                  std::list<BCError>&      theErrors)
2127 {
2128   SetErrorCode(KO);
2129
2130   if (theCompound.IsNull()) return Standard_False;
2131   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2132
2133   Standard_Boolean isCompOfBlocks = Standard_True;
2134
2135   // Map sub-shapes and their indices
2136   TopTools_IndexedMapOfShape anIndices;
2137   TopExp::MapShapes(aBlockOrComp, anIndices);
2138
2139   // 1. Report non-blocks
2140   TopTools_ListOfShape NOT; // Not blocks
2141   TopTools_ListOfShape DEG; // Hexahedral solids, having degenerated edges
2142   TopTools_ListOfShape SEA; // Hexahedral solids, having seam edges
2143   TopTools_ListOfShape BLO; // All blocks from the given compound
2144   AddBlocksFromOld(aBlockOrComp, BLO, NOT, DEG, SEA);
2145
2146   if (NOT.Extent() > 0) {
2147     isCompOfBlocks = Standard_False;
2148     BCError anErr;
2149     anErr.error = NOT_BLOCK;
2150     TopTools_ListIteratorOfListOfShape it (NOT);
2151     for (; it.More(); it.Next()) {
2152       anErr.incriminated.push_back(anIndices.FindIndex(it.Value()));
2153     }
2154     theErrors.push_back(anErr);
2155   }
2156
2157   if (DEG.Extent() > 0 || SEA.Extent() > 0) {
2158     isCompOfBlocks = Standard_False;
2159     BCError anErr;
2160     anErr.error = EXTRA_EDGE;
2161
2162     TopTools_ListIteratorOfListOfShape itDEG (DEG);
2163     for (; itDEG.More(); itDEG.Next()) {
2164       anErr.incriminated.push_back(anIndices.FindIndex(itDEG.Value()));
2165     }
2166
2167     TopTools_ListIteratorOfListOfShape itSEA (SEA);
2168     for (; itSEA.More(); itSEA.Next()) {
2169       anErr.incriminated.push_back(anIndices.FindIndex(itSEA.Value()));
2170     }
2171
2172     theErrors.push_back(anErr);
2173   }
2174
2175   Standard_Integer nbBlocks = BLO.Extent();
2176   if (nbBlocks == 0) {
2177     isCompOfBlocks = Standard_False;
2178     SetErrorCode(OK);
2179     return isCompOfBlocks;
2180   }
2181   if (nbBlocks == 1) {
2182     SetErrorCode(OK);
2183     return isCompOfBlocks;
2184   }
2185
2186   // Convert list of blocks into array for easy and fast access
2187   Standard_Integer ibl = 1;
2188   TopTools_Array1OfShape aBlocks (1, nbBlocks);
2189   TopTools_ListIteratorOfListOfShape BLOit (BLO);
2190   for (; BLOit.More(); BLOit.Next(), ibl++) {
2191     aBlocks.SetValue(ibl, BLOit.Value());
2192   }
2193
2194   // 2. Find relations between all blocks,
2195   //    report connection errors (NOT_GLUED and INVALID_CONNECTION)
2196   TColStd_Array2OfInteger aRelations (1, nbBlocks, 1, nbBlocks);
2197   aRelations.Init(REL_NOT_CONNECTED);
2198
2199   Standard_Integer row = 1;
2200   for (row = 1; row <= nbBlocks; row++) {
2201     TopoDS_Shape aBlock = aBlocks.Value(row);
2202
2203     Standard_Integer col = row + 1;
2204     for (; col <= nbBlocks; col++) {
2205       Standard_Integer aRel = BlocksRelation(aBlock, aBlocks.Value(col));
2206       if (aRel != REL_NOT_CONNECTED) {
2207         aRelations.SetValue(row, col, aRel);
2208         aRelations.SetValue(col, row, aRel);
2209         if (aRel == REL_NOT_GLUED) {
2210           // report connection error
2211           isCompOfBlocks = Standard_False;
2212           BCError anErr;
2213           anErr.error = NOT_GLUED;
2214           anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(row)));
2215           anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(col)));
2216           theErrors.push_back(anErr);
2217         } else if (aRel == REL_COLLISION_VV ||
2218                    aRel == REL_COLLISION_FF ||
2219                    aRel == REL_COLLISION_EE ||
2220                    aRel == REL_UNKNOWN) {
2221           // report connection error
2222           isCompOfBlocks = Standard_False;
2223           BCError anErr;
2224           anErr.error = INVALID_CONNECTION;
2225           anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(row)));
2226           anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(col)));
2227           theErrors.push_back(anErr);
2228         } else {
2229         }
2230       }
2231     }
2232   }
2233
2234   // 3. Find largest set of connected (good connection or not glued) blocks
2235   TColStd_MapOfInteger aProcessedMap;
2236   TColStd_MapOfInteger aLargestSet;
2237   TColStd_MapOfInteger aCurrentSet;
2238   for (ibl = 1; ibl <= nbBlocks; ibl++) {
2239     if (!aProcessedMap.Contains(ibl)) {
2240       aCurrentSet.Clear();
2241       FindConnected(ibl, aRelations, aProcessedMap, aCurrentSet);
2242       if (aCurrentSet.Extent() > aLargestSet.Extent()) {
2243         aLargestSet = aCurrentSet;
2244       }
2245     }
2246   }
2247
2248   // 4. Report all blocks, isolated from <aLargestSet>
2249   BCError anErr;
2250   anErr.error = NOT_CONNECTED;
2251   Standard_Boolean hasIsolated = Standard_False;
2252   for (ibl = 1; ibl <= nbBlocks; ibl++) {
2253     if (!aLargestSet.Contains(ibl)) {
2254       aProcessedMap.Clear();
2255       if (!HasAnyConnection(ibl, aLargestSet, aRelations, aProcessedMap)) {
2256         // report connection absence
2257         hasIsolated = Standard_True;
2258         anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(ibl)));
2259       }
2260     }
2261   }
2262   if (hasIsolated) {
2263     isCompOfBlocks = Standard_False;
2264     theErrors.push_back(anErr);
2265   }
2266
2267   SetErrorCode(OK);
2268   return isCompOfBlocks;
2269 }
2270
2271 //=============================================================================
2272 /*!
2273  *  PrintBCErrors
2274  */
2275 //=============================================================================
2276 TCollection_AsciiString GEOMImpl_IBlocksOperations::PrintBCErrors
2277                                                 (Handle(GEOM_Object)  theCompound,
2278                                                  const std::list<BCError>& theErrors)
2279 {
2280   TCollection_AsciiString aDescr;
2281
2282   std::list<BCError>::const_iterator errIt = theErrors.begin();
2283   int i = 0;
2284   for (; errIt != theErrors.end(); i++, errIt++) {
2285     BCError errStruct = *errIt;
2286
2287     switch (errStruct.error) {
2288     case NOT_BLOCK:
2289       aDescr += "\n\tNot a Blocks: ";
2290       break;
2291     case EXTRA_EDGE:
2292       aDescr += "\n\tHexahedral solids with degenerated and/or seam edges: ";
2293       break;
2294     case INVALID_CONNECTION:
2295       aDescr += "\n\tInvalid connection between two blocks: ";
2296       break;
2297     case NOT_CONNECTED:
2298       aDescr += "\n\tBlocks, not connected with main body: ";
2299       break;
2300     case NOT_GLUED:
2301       aDescr += "\n\tNot glued blocks: ";
2302       break;
2303     default:
2304       break;
2305     }
2306
2307     std::list<int> sshList = errStruct.incriminated;
2308     std::list<int>::iterator sshIt = sshList.begin();
2309     int jj = 0;
2310     for (; sshIt != sshList.end(); jj++, sshIt++) {
2311       if (jj > 0)
2312         aDescr += ", ";
2313       aDescr += TCollection_AsciiString(*sshIt);
2314     }
2315   }
2316
2317   return aDescr;
2318 }
2319
2320 //=============================================================================
2321 /*!
2322  *  CheckCompoundOfBlocks
2323  */
2324 //=============================================================================
2325 Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
2326                                               (Handle(GEOM_Object) theCompound,
2327                                                std::list<BCError>& theErrors)
2328 {
2329   SetErrorCode(KO);
2330
2331   if (theCompound.IsNull()) return Standard_False;
2332   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2333
2334   Standard_Boolean isCompOfBlocks = Standard_True;
2335
2336   // Map sub-shapes and their indices
2337   TopTools_IndexedMapOfShape anIndices;
2338   TopExp::MapShapes(aBlockOrComp, anIndices);
2339
2340   // 1. Separate blocks from non-blocks
2341   TopTools_ListOfShape NOT; // Not blocks
2342   TopTools_ListOfShape EXT; // Hexahedral solids, having degenerated and/or seam edges
2343   TopTools_ListOfShape BLO; // All blocks from the given compound
2344   TopTools_ListOfShape NOQ; // All non-quadrangular faces
2345   AddBlocksFrom(aBlockOrComp, BLO, NOT, EXT, NOQ);
2346
2347   // Report non-blocks
2348   if (NOT.Extent() > 0) {
2349     isCompOfBlocks = Standard_False;
2350     BCError anErr;
2351     anErr.error = NOT_BLOCK;
2352     TopTools_ListIteratorOfListOfShape it (NOT);
2353     for (; it.More(); it.Next()) {
2354       anErr.incriminated.push_back(anIndices.FindIndex(it.Value()));
2355     }
2356     theErrors.push_back(anErr);
2357   }
2358
2359   // Report solids, having degenerated and/or seam edges
2360   if (EXT.Extent() > 0) {
2361     isCompOfBlocks = Standard_False;
2362     BCError anErr;
2363     anErr.error = EXTRA_EDGE;
2364     TopTools_ListIteratorOfListOfShape it (EXT);
2365     for (; it.More(); it.Next()) {
2366       anErr.incriminated.push_back(anIndices.FindIndex(it.Value()));
2367     }
2368     theErrors.push_back(anErr);
2369   }
2370
2371   Standard_Integer nbBlocks = BLO.Extent();
2372   if (nbBlocks == 0) {
2373     isCompOfBlocks = Standard_False;
2374     SetErrorCode(OK);
2375     return isCompOfBlocks;
2376   }
2377   if (nbBlocks == 1) {
2378     SetErrorCode(OK);
2379     return isCompOfBlocks;
2380   }
2381
2382   // Prepare data for 2. and 3.
2383   TColStd_Array2OfInteger aRelations (1, nbBlocks, 1, nbBlocks);
2384   aRelations.Init(REL_NOT_CONNECTED);
2385
2386   TopTools_IndexedMapOfShape mapBlocks;
2387
2388   BRep_Builder BB;
2389   TopoDS_Compound aComp;
2390   BB.MakeCompound(aComp);
2391
2392   TopTools_ListIteratorOfListOfShape BLOit (BLO);
2393   for (; BLOit.More(); BLOit.Next()) {
2394     mapBlocks.Add(BLOit.Value());
2395     BB.Add(aComp, BLOit.Value());
2396   }
2397
2398   // 2. Find glued blocks (having shared faces)
2399   TopTools_IndexedDataMapOfShapeListOfShape mapFaceBlocks;
2400   GEOMImpl_Block6Explorer::MapShapesAndAncestors
2401     (aComp, TopAbs_FACE, TopAbs_SOLID, mapFaceBlocks);
2402
2403   Standard_Integer prevInd = 0, curInd = 0;
2404   Standard_Integer ind = 1, nbFaces = mapFaceBlocks.Extent();
2405   for (; ind <= nbFaces; ind++) {
2406     const TopTools_ListOfShape& aGluedBlocks = mapFaceBlocks.FindFromIndex(ind);
2407     if (aGluedBlocks.Extent() > 1) { // Shared face found
2408       TopTools_ListIteratorOfListOfShape aGluedBlocksIt (aGluedBlocks);
2409       TopoDS_Shape prevBlock, curBlock;
2410       for (; aGluedBlocksIt.More(); aGluedBlocksIt.Next()) {
2411         curBlock = aGluedBlocksIt.Value();
2412         if (!prevBlock.IsNull()) {
2413           prevInd = mapBlocks.FindIndex(prevBlock);
2414           curInd  = mapBlocks.FindIndex(curBlock);
2415           aRelations.SetValue(prevInd, curInd, REL_OK);
2416           aRelations.SetValue(curInd, prevInd, REL_OK);
2417         }
2418         prevBlock = curBlock;
2419       }
2420     }
2421   }
2422
2423   // 3. Find not glued blocks
2424   GEOMAlgo_GlueAnalyser aGD;
2425
2426   aGD.SetShape(aComp);
2427   aGD.SetTolerance(Precision::Confusion());
2428   aGD.SetCheckGeometry(Standard_True);
2429   aGD.Perform();
2430
2431   Standard_Integer iErr, iWrn;
2432   iErr = aGD.ErrorStatus();
2433   if (iErr) {
2434     SetErrorCode("Error in GEOMAlgo_GlueAnalyser");
2435     return isCompOfBlocks;
2436   }
2437   iWrn = aGD.WarningStatus();
2438   if (iWrn) {
2439     MESSAGE("Warning in GEOMAlgo_GlueAnalyser");
2440   }
2441
2442   // Report not glued blocks
2443   if (aGD.HasSolidsToGlue()) {
2444     isCompOfBlocks = Standard_False;
2445     Standard_Integer aSx1Ind, aSx2Ind;
2446
2447     const GEOMAlgo_ListOfCoupleOfShapes& aLCS = aGD.SolidsToGlue();
2448     GEOMAlgo_ListIteratorOfListOfCoupleOfShapes aItCS (aLCS);
2449     for (; aItCS.More(); aItCS.Next()) {
2450       const GEOMAlgo_CoupleOfShapes& aCS = aItCS.Value();
2451       const TopoDS_Shape& aSx1 = aCS.Shape1();
2452       const TopoDS_Shape& aSx2 = aCS.Shape2();
2453
2454       aSx1Ind = mapBlocks.FindIndex(aSx1);
2455       aSx2Ind = mapBlocks.FindIndex(aSx2);
2456       aRelations.SetValue(aSx1Ind, aSx2Ind, NOT_GLUED);
2457       aRelations.SetValue(aSx2Ind, aSx1Ind, NOT_GLUED);
2458
2459       BCError anErr;
2460       anErr.error = NOT_GLUED;
2461       anErr.incriminated.push_back(anIndices.FindIndex(aSx1));
2462       anErr.incriminated.push_back(anIndices.FindIndex(aSx2));
2463       theErrors.push_back(anErr);
2464     }
2465   }
2466
2467   // 4. Find largest set of connected (good connection or not glued) blocks
2468   Standard_Integer ibl = 1;
2469   TColStd_MapOfInteger aProcessedMap;
2470   TColStd_MapOfInteger aLargestSet;
2471   TColStd_MapOfInteger aCurrentSet;
2472   for (ibl = 1; ibl <= nbBlocks; ibl++) {
2473     if (!aProcessedMap.Contains(ibl)) {
2474       aCurrentSet.Clear();
2475       FindConnected(ibl, aRelations, aProcessedMap, aCurrentSet);
2476       if (aCurrentSet.Extent() > aLargestSet.Extent()) {
2477         aLargestSet = aCurrentSet;
2478       }
2479     }
2480   }
2481
2482   // 5. Report all blocks, isolated from <aLargestSet>
2483   BCError anErr;
2484   anErr.error = NOT_CONNECTED;
2485   Standard_Boolean hasIsolated = Standard_False;
2486   for (ibl = 1; ibl <= nbBlocks; ibl++) {
2487     if (!aLargestSet.Contains(ibl)) {
2488       aProcessedMap.Clear();
2489       if (!HasAnyConnection(ibl, aLargestSet, aRelations, aProcessedMap)) {
2490         // report connection absence
2491         hasIsolated = Standard_True;
2492         anErr.incriminated.push_back(anIndices.FindIndex(mapBlocks.FindKey(ibl)));
2493       }
2494     }
2495   }
2496   if (hasIsolated) {
2497     isCompOfBlocks = Standard_False;
2498     theErrors.push_back(anErr);
2499   }
2500
2501   SetErrorCode(OK);
2502   return isCompOfBlocks;
2503 }
2504
2505 //=============================================================================
2506 /*!
2507  *  GetNonBlocks
2508  */
2509 //=============================================================================
2510 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetNonBlocks
2511                                      (Handle(GEOM_Object) theShape,
2512                                       Handle(GEOM_Object)& theNonQuads)
2513 {
2514   SetErrorCode(KO);
2515
2516   if (theShape.IsNull()) return NULL;
2517   TopoDS_Shape aShape = theShape->GetValue();
2518
2519   // Separate blocks from non-blocks
2520   TopTools_ListOfShape BLO; // All blocks from the given compound
2521   TopTools_ListOfShape NOT; // Not blocks
2522   TopTools_ListOfShape EXT; // Hexahedral solids, having degenerated and/or seam edges
2523   TopTools_ListOfShape NOQ; // All non-quadrangular faces
2524   AddBlocksFrom(aShape, BLO, NOT, EXT, NOQ);
2525
2526   if (NOT.IsEmpty() && EXT.IsEmpty() && NOQ.IsEmpty()) {
2527     SetErrorCode("NOT_FOUND_ANY");
2528     return NULL;
2529   }
2530
2531   // Map sub-shapes and their indices
2532   TopTools_IndexedMapOfShape anIndices;
2533   TopExp::MapShapes(aShape, anIndices);
2534
2535   // Non-blocks
2536   Handle(GEOM_Object) aNonBlocks;
2537   if (NOT.Extent() > 0 || EXT.Extent() > 0) {
2538     Handle(TColStd_HArray1OfInteger) anArray =
2539       new TColStd_HArray1OfInteger (1, NOT.Extent() + EXT.Extent());
2540     Standard_Integer ii = 1;
2541     TopTools_ListIteratorOfListOfShape it1 (NOT);
2542     for (; it1.More(); it1.Next(), ii++) {
2543       anArray->SetValue(ii, anIndices.FindIndex(it1.Value()));
2544     }
2545     TopTools_ListIteratorOfListOfShape it2 (EXT);
2546     for (; it2.More(); it2.Next(), ii++) {
2547       anArray->SetValue(ii, anIndices.FindIndex(it2.Value()));
2548     }
2549
2550     aNonBlocks = GetEngine()->AddSubShape(theShape, anArray);
2551     if (aNonBlocks.IsNull()) {
2552       SetErrorCode("Error in algorithm: result found, but cannot be returned.");
2553       return NULL;
2554     }
2555     aNonBlocks->SetType(GEOM_GROUP);
2556     TDF_Label aFreeLabel = aNonBlocks->GetFreeLabel();
2557     TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)TopAbs_SOLID);
2558   }
2559
2560   // Non-quadrangles
2561   if (NOQ.Extent() > 0) {
2562     Handle(TColStd_HArray1OfInteger) anArray =
2563       new TColStd_HArray1OfInteger (1, NOQ.Extent());
2564     Standard_Integer ii = 1;
2565     TopTools_ListIteratorOfListOfShape it1 (NOQ);
2566     for (; it1.More(); it1.Next(), ii++) {
2567       anArray->SetValue(ii, anIndices.FindIndex(it1.Value()));
2568     }
2569
2570     theNonQuads = GetEngine()->AddSubShape(theShape, anArray);
2571     if (theNonQuads.IsNull()) {
2572       SetErrorCode("Error in algorithm: result found, but cannot be returned.");
2573       return NULL;
2574     }
2575     theNonQuads->SetType(GEOM_GROUP);
2576     TDF_Label aFreeLabel = theNonQuads->GetFreeLabel();
2577     TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)TopAbs_FACE);
2578   }
2579
2580   //Make a Python command
2581   Handle(GEOM_Function) aMainShape = theShape->GetLastFunction();
2582   GEOM::TPythonDump pd (aMainShape, /*append=*/true);
2583   pd << "(";
2584   if (aNonBlocks.IsNull())
2585     pd << "no_bad_solids";
2586   else
2587     pd << aNonBlocks;
2588   pd << ", ";
2589   if (theNonQuads.IsNull())
2590     pd << "no_bad_faces";
2591   else
2592     pd << theNonQuads;
2593   pd << ") = geompy.GetNonBlocks(" << theShape << ")";
2594
2595   SetErrorCode(OK);
2596   return aNonBlocks;
2597 }
2598
2599 //=============================================================================
2600 /*!
2601  *  RemoveExtraEdges
2602  */
2603 //=============================================================================
2604 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::RemoveExtraEdges
2605                                      (Handle(GEOM_Object) theObject,
2606                                       const Standard_Integer theOptimumNbFaces)
2607 {
2608   SetErrorCode(KO);
2609
2610   if (theObject.IsNull()) return NULL;
2611
2612   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
2613   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be fixed
2614
2615   //Add a new Copy object
2616   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
2617
2618   //Add a function
2619   Handle(GEOM_Function) aFunction =
2620     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_REMOVE_EXTRA);
2621
2622   //Check if the function is set correctly
2623   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
2624
2625   GEOMImpl_IBlockTrsf aTI (aFunction);
2626   aTI.SetOriginal(aLastFunction);
2627   aTI.SetOptimumNbFaces(theOptimumNbFaces);
2628
2629   //Compute the fixed shape
2630   try {
2631 #if OCC_VERSION_LARGE > 0x06010000
2632     OCC_CATCH_SIGNALS;
2633 #endif
2634     if (!GetSolver()->ComputeFunction(aFunction)) {
2635       SetErrorCode("Block driver failed to remove extra edges of the given shape");
2636       return NULL;
2637     }
2638   }
2639   catch (Standard_Failure) {
2640     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2641     SetErrorCode(aFail->GetMessageString());
2642     return NULL;
2643   }
2644
2645   //Make a Python command
2646   std::string doUnionFaces = (theOptimumNbFaces < 0) ? "False" : "True";
2647   GEOM::TPythonDump(aFunction) << aCopy << " = geompy.RemoveExtraEdges("
2648                                << theObject << ", " << doUnionFaces.data() << ")";
2649
2650   SetErrorCode(OK);
2651   return aCopy;
2652 }
2653
2654 //=============================================================================
2655 /*!
2656  *  UnionFaces
2657  */
2658 //=============================================================================
2659 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::UnionFaces
2660                                      (Handle(GEOM_Object) theObject)
2661 {
2662   SetErrorCode(KO);
2663
2664   if (theObject.IsNull()) return NULL;
2665
2666   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
2667   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be fixed
2668
2669   //Add a new Copy object
2670   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
2671
2672   //Add a function
2673   Handle(GEOM_Function) aFunction =
2674     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_UNION_FACES);
2675
2676   //Check if the function is set correctly
2677   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
2678
2679   GEOMImpl_IBlockTrsf aTI (aFunction);
2680   aTI.SetOriginal(aLastFunction);
2681
2682   //Compute the fixed shape
2683   try {
2684 #if OCC_VERSION_LARGE > 0x06010000
2685     OCC_CATCH_SIGNALS;
2686 #endif
2687     if (!GetSolver()->ComputeFunction(aFunction)) {
2688       SetErrorCode("Block driver failed to remove extra edges of the given shape");
2689       return NULL;
2690     }
2691   }
2692   catch (Standard_Failure) {
2693     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2694     SetErrorCode(aFail->GetMessageString());
2695     return NULL;
2696   }
2697
2698   //Make a Python command
2699   GEOM::TPythonDump(aFunction) << aCopy << " = geompy.UnionFaces("
2700                                << theObject << ")";
2701
2702   SetErrorCode(OK);
2703   return aCopy;
2704 }
2705
2706 //=============================================================================
2707 /*!
2708  *  CheckAndImprove
2709  */
2710 //=============================================================================
2711 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::CheckAndImprove
2712                                              (Handle(GEOM_Object) theObject)
2713 {
2714   SetErrorCode(KO);
2715
2716   if (theObject.IsNull()) return NULL;
2717
2718   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
2719   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be fixed
2720
2721   //Add a new Copy object
2722   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
2723
2724   //Add a function
2725   Handle(GEOM_Function) aFunction =
2726     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_COMPOUND_IMPROVE);
2727
2728   //Check if the function is set correctly
2729   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
2730
2731   GEOMImpl_IBlockTrsf aTI (aFunction);
2732   aTI.SetOriginal(aLastFunction);
2733
2734   // -1 means do not unite faces on common surface (?except case of seam edge between them?)
2735   //aTI.SetOptimumNbFaces(-1);
2736   aTI.SetOptimumNbFaces(6);
2737
2738   //Compute the fixed shape
2739   try {
2740 #if OCC_VERSION_LARGE > 0x06010000
2741     OCC_CATCH_SIGNALS;
2742 #endif
2743     if (!GetSolver()->ComputeFunction(aFunction)) {
2744       SetErrorCode("Block driver failed to improve the given blocks compound");
2745       return NULL;
2746     }
2747   }
2748   catch (Standard_Failure) {
2749     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2750     SetErrorCode(aFail->GetMessageString());
2751     return NULL;
2752   }
2753
2754   //Make a Python command
2755   GEOM::TPythonDump(aFunction) << aCopy
2756     << " = geompy.CheckAndImprove(" << theObject << ")";
2757
2758   SetErrorCode(OK);
2759   return aCopy;
2760 }
2761
2762 //=============================================================================
2763 /*!
2764  *  ExplodeCompoundOfBlocks
2765  */
2766 //=============================================================================
2767 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::ExplodeCompoundOfBlocks
2768                                                 (Handle(GEOM_Object)    theCompound,
2769                                                  const Standard_Integer theMinNbFaces,
2770                                                  const Standard_Integer theMaxNbFaces)
2771 {
2772   SetErrorCode(KO);
2773
2774   if (theCompound.IsNull()) return NULL;
2775   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2776   if (aBlockOrComp.IsNull()) return NULL;
2777
2778   Handle(TColStd_HSequenceOfTransient) aBlocks = new TColStd_HSequenceOfTransient;
2779   Handle(GEOM_Object) anObj;
2780   Handle(GEOM_Function) aFunction;
2781
2782   TopTools_MapOfShape mapShape;
2783   TCollection_AsciiString anAsciiList, anEntry;
2784
2785   // Map shapes
2786   TopTools_IndexedMapOfShape anIndices;
2787   TopExp::MapShapes(aBlockOrComp, anIndices);
2788   Handle(TColStd_HArray1OfInteger) anArray;
2789
2790   // Explode
2791   try {
2792 #if OCC_VERSION_LARGE > 0x06010000
2793     OCC_CATCH_SIGNALS;
2794 #endif
2795     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
2796     for (; exp.More(); exp.Next()) {
2797       if (mapShape.Add(exp.Current())) {
2798         TopoDS_Shape aSolid = exp.Current();
2799
2800         TopTools_MapOfShape mapFaces;
2801         TopExp_Explorer expF (aSolid, TopAbs_FACE);
2802         Standard_Integer nbFaces = 0;
2803         for (; expF.More(); expF.Next()) {
2804           if (mapFaces.Add(expF.Current())) {
2805             nbFaces++;
2806           }
2807         }
2808
2809         if (theMinNbFaces <= nbFaces && nbFaces <= theMaxNbFaces) {
2810           anArray = new TColStd_HArray1OfInteger(1,1);
2811           anArray->SetValue(1, anIndices.FindIndex(aSolid));
2812           anObj = GetEngine()->AddSubShape(theCompound, anArray);
2813           aBlocks->Append(anObj);
2814
2815           //Make a Python command
2816           TDF_Tool::Entry(anObj->GetEntry(), anEntry);
2817           anAsciiList += anEntry + ", ";
2818         }
2819       }
2820     }
2821   }
2822   catch (Standard_Failure) {
2823     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2824     SetErrorCode(aFail->GetMessageString());
2825     return aBlocks;
2826   }
2827
2828   if (aBlocks->IsEmpty()) {
2829     SetErrorCode("There are no specified blocks in the given shape");
2830     return aBlocks;
2831   }
2832
2833   anAsciiList.Trunc(anAsciiList.Length() - 2);
2834
2835   //The explode doesn't change object so no new function is required.
2836   aFunction = theCompound->GetLastFunction();
2837
2838   //Make a Python command
2839   GEOM::TPythonDump(aFunction, /*append=*/true)
2840     << "[" << anAsciiList.ToCString() << "] = geompy.MakeBlockExplode("
2841     << theCompound << ", " << theMinNbFaces << ", " << theMaxNbFaces << ")";
2842
2843   SetErrorCode(OK);
2844   return aBlocks;
2845 }
2846
2847 //=============================================================================
2848 /*!
2849  *  GetBlockNearPoint
2850  */
2851 //=============================================================================
2852 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockNearPoint
2853                                                 (Handle(GEOM_Object) theCompound,
2854                                                  Handle(GEOM_Object) thePoint)
2855 {
2856   SetErrorCode(KO);
2857
2858   //New object
2859   Handle(GEOM_Object) aResult;
2860
2861   // Arguments
2862   if (theCompound.IsNull() || thePoint.IsNull()) return NULL;
2863
2864   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2865   if (aBlockOrComp.IsNull()) {
2866     SetErrorCode("Compound is null");
2867     return NULL;
2868   }
2869   if (aBlockOrComp.ShapeType() != TopAbs_COMPOUND &&
2870       aBlockOrComp.ShapeType() != TopAbs_COMPSOLID) {
2871     SetErrorCode("Shape to find block in is not a compound");
2872     return NULL;
2873   }
2874
2875   TopoDS_Shape anArg = thePoint->GetValue();
2876   if (anArg.IsNull()) {
2877     SetErrorCode("Point is null");
2878     return NULL;
2879   }
2880   if (anArg.ShapeType() != TopAbs_VERTEX) {
2881     SetErrorCode("Shape for block identification is not a vertex");
2882     return NULL;
2883   }
2884
2885   //Compute the Block value
2886   try {
2887 #if OCC_VERSION_LARGE > 0x06010000
2888     OCC_CATCH_SIGNALS;
2889 #endif
2890     TopoDS_Shape aShape;
2891
2892     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
2893     gp_Pnt aPnt = BRep_Tool::Pnt(aVert);
2894     Standard_Real PX, PY, PZ;
2895     aPnt.Coord(PX, PY, PZ);
2896
2897     // 1. Classify the point relatively each block
2898     Standard_Integer nearest = 2, nbFound = 0;
2899     TopTools_DataMapOfShapeInteger mapShapeDist;
2900     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
2901     for (; exp.More(); exp.Next()) {
2902       TopoDS_Shape aSolid = exp.Current();
2903
2904       if (!mapShapeDist.IsBound(aSolid)) {
2905         Standard_Integer aDistance = 2;
2906
2907         // 1.a. Classify relatively Bounding box
2908         Standard_Real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
2909         Bnd_Box BB;
2910         BRepBndLib::Add(aSolid, BB);
2911         BB.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
2912         if (PX < Xmin || Xmax < PX ||
2913             PY < Ymin || Ymax < PY ||
2914             PZ < Zmin || Zmax < PZ) {
2915           // OUT of bounding box
2916           aDistance = 1;
2917         } else {
2918           // 1.b. Classify relatively the solid itself
2919           BRepClass3d_SolidClassifier SC (aSolid, aPnt, Precision::Confusion());
2920           if (SC.State() == TopAbs_IN) {
2921             aDistance = -1;
2922           } else if (SC.State() == TopAbs_ON) {
2923             aDistance = 0;
2924           } else { // OUT
2925             aDistance = 1;
2926           }
2927         }
2928
2929         if (aDistance < nearest) {
2930           nearest = aDistance;
2931           aShape = aSolid;
2932           nbFound = 1;
2933
2934           // A first found block, containing the point inside, will be returned.
2935           // It is the solution, if there are no intersecting blocks in the compound.
2936           if (nearest == -1) break;
2937
2938         } else if (aDistance == nearest) {
2939           nbFound++;
2940         } else {
2941         }
2942
2943         mapShapeDist.Bind(aSolid, aDistance);
2944       } // if (!mapShapeDist.IsBound(aSolid))
2945     }
2946
2947     // 2. Define block, containing the point or having minimum distance to it
2948     if (nbFound > 1) {
2949       if (nearest == 0) {
2950         // The point is on boundary of some blocks and there are
2951         // no blocks, having the point inside their volume
2952         SetErrorCode("Multiple blocks near the given point are found");
2953         return NULL;
2954
2955       } else if (nearest == 1) {
2956         // The point is outside some blocks and there are
2957         // no blocks, having the point inside or on boundary.
2958         // We will get a nearest block
2959         Standard_Real minDist = RealLast();
2960         TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
2961         for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
2962           if (mapShapeDistIter.Value() == 1) {
2963             TopoDS_Shape aSolid = mapShapeDistIter.Key();
2964             BRepExtrema_DistShapeShape aDistTool (aVert, aSolid);
2965             if (!aDistTool.IsDone()) {
2966               SetErrorCode("Can not find a distance from the given point to one of blocks");
2967               return NULL;
2968             }
2969             Standard_Real aDist = aDistTool.Value();
2970             if (aDist < minDist) {
2971               minDist = aDist;
2972               aShape = aSolid;
2973             }
2974           }
2975         }
2976       } else { // nearest == -1
2977 //        // The point is inside some blocks.
2978 //        // We will get a block with nearest center
2979 //        Standard_Real minDist = RealLast();
2980 //        TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
2981 //        for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
2982 //          if (mapShapeDistIter.Value() == -1) {
2983 //            TopoDS_Shape aSolid = mapShapeDistIter.Key();
2984 //            GProp_GProps aSystem;
2985 //            BRepGProp::VolumeProperties(aSolid, aSystem);
2986 //            gp_Pnt aCenterMass = aSystem.CentreOfMass();
2987 //
2988 //            Standard_Real aDist = aCenterMass.Distance(aPnt);
2989 //            if (aDist < minDist) {
2990 //              minDist = aDist;
2991 //              aShape = aSolid;
2992 //            }
2993 //          }
2994 //        }
2995       }
2996     } // if (nbFound > 1)
2997
2998     if (nbFound == 0) {
2999       SetErrorCode("There are no blocks near the given point");
3000       return NULL;
3001     } else {
3002       TopTools_IndexedMapOfShape anIndices;
3003       TopExp::MapShapes(aBlockOrComp, anIndices);
3004       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
3005       anArray->SetValue(1, anIndices.FindIndex(aShape));
3006       aResult = GetEngine()->AddSubShape(theCompound, anArray);
3007     }
3008   }
3009   catch (Standard_Failure) {
3010     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3011     SetErrorCode(aFail->GetMessageString());
3012     return NULL;
3013   }
3014
3015   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
3016
3017   //Make a Python command
3018   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetBlockNearPoint("
3019                                << theCompound << ", " << thePoint << ")";
3020
3021   SetErrorCode(OK);
3022   return aResult;
3023 }
3024
3025 //=============================================================================
3026 /*!
3027  *  GetBlockByParts
3028  */
3029 //=============================================================================
3030 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockByParts
3031                       (Handle(GEOM_Object)                         theCompound,
3032                        const Handle(TColStd_HSequenceOfTransient)& theParts)
3033 {
3034   SetErrorCode(KO);
3035
3036   Handle(GEOM_Object) aResult;
3037
3038   if (theCompound.IsNull() || theParts.IsNull()) return NULL;
3039   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
3040   if (aBlockOrComp.IsNull()) return NULL;
3041
3042   //Get the parts
3043   Standard_Integer argi, aLen = theParts->Length();
3044   TopTools_Array1OfShape anArgs (1, aLen);
3045   TCollection_AsciiString anEntry, aPartsDescr;
3046   for (argi = 1; argi <= aLen; argi++) {
3047     Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast(theParts->Value(argi));
3048     Handle(GEOM_Function) aRef = anObj->GetLastFunction();
3049     if (aRef.IsNull()) return NULL;
3050
3051     TopoDS_Shape anArg = aRef->GetValue();
3052     if (anArg.IsNull()) {
3053       SetErrorCode("Null shape is given as argument");
3054       return NULL;
3055     }
3056     anArgs(argi) = anArg;
3057
3058     // For Python command
3059     TDF_Tool::Entry(anObj->GetEntry(), anEntry);
3060     if (argi > 1) aPartsDescr += ", ";
3061     aPartsDescr += anEntry;
3062   }
3063
3064   //Compute the Block value
3065   try {
3066 #if OCC_VERSION_LARGE > 0x06010000
3067     OCC_CATCH_SIGNALS;
3068 #endif
3069     // 1. Explode compound on solids
3070     TopTools_MapOfShape mapShape;
3071     Standard_Integer nbSolids = 0;
3072     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
3073     for (; exp.More(); exp.Next()) {
3074       if (mapShape.Add(exp.Current())) {
3075         nbSolids++;
3076       }
3077     }
3078
3079     mapShape.Clear();
3080     Standard_Integer ind = 1;
3081     TopTools_Array1OfShape aSolids (1, nbSolids);
3082     TColStd_Array1OfInteger aNbParts (1, nbSolids);
3083     for (exp.Init(aBlockOrComp, TopAbs_SOLID); exp.More(); exp.Next(), ind++) {
3084       if (mapShape.Add(exp.Current())) {
3085         TopoDS_Shape aSolid = exp.Current();
3086         aSolids(ind) = aSolid;
3087         aNbParts(ind) = 0;
3088
3089         // 2. Define quantity of parts, contained in each solid
3090         TopTools_IndexedMapOfShape aSubShapes;
3091         TopExp::MapShapes(aSolid, aSubShapes);
3092         for (argi = 1; argi <= aLen; argi++) {
3093           if (aSubShapes.Contains(anArgs(argi))) {
3094             aNbParts(ind)++;
3095           }
3096         }
3097       }
3098     }
3099
3100     // 3. Define solid, containing maximum quantity of parts
3101     Standard_Integer maxNb = 0, nbFound = 0;
3102     TopoDS_Shape aShape;
3103     for (ind = 1; ind <= nbSolids; ind++) {
3104       if (aNbParts(ind) > maxNb) {
3105         maxNb = aNbParts(ind);
3106         aShape = aSolids(ind);
3107         nbFound = 1;
3108       } else if (aNbParts(ind) == maxNb) {
3109         nbFound++;
3110       } else {
3111       }
3112     }
3113     if (nbFound > 1) {
3114       SetErrorCode("Multiple blocks, containing maximum quantity of the given parts, are found");
3115       return NULL;
3116     } else if (nbFound == 0) {
3117       SetErrorCode("There are no blocks, containing the given parts");
3118       return NULL;
3119     } else {
3120       TopTools_IndexedMapOfShape anIndices;
3121       TopExp::MapShapes(aBlockOrComp, anIndices);
3122       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
3123       anArray->SetValue(1, anIndices.FindIndex(aShape));
3124       aResult = GetEngine()->AddSubShape(theCompound, anArray);
3125     }
3126   } catch (Standard_Failure) {
3127     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3128     SetErrorCode(aFail->GetMessageString());
3129     return NULL;
3130   }
3131
3132   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
3133
3134   //Make a Python command
3135   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetBlockByParts("
3136     << theCompound << ", [" << aPartsDescr.ToCString() << "])";
3137
3138   SetErrorCode(OK);
3139   return aResult;
3140 }
3141
3142 //=============================================================================
3143 /*!
3144  *  GetBlocksByParts
3145  */
3146 //=============================================================================
3147 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::GetBlocksByParts
3148                       (Handle(GEOM_Object)                         theCompound,
3149                        const Handle(TColStd_HSequenceOfTransient)& theParts)
3150 {
3151   SetErrorCode(KO);
3152
3153   if (theCompound.IsNull() || theParts.IsNull()) return NULL;
3154   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
3155   if (aBlockOrComp.IsNull()) return NULL;
3156
3157   Handle(TColStd_HSequenceOfTransient) aBlocks = new TColStd_HSequenceOfTransient;
3158   Handle(GEOM_Object) anObj;
3159   Handle(GEOM_Function) aFunction;
3160
3161   //Get the parts
3162   Standard_Integer argi, aLen = theParts->Length();
3163   TopTools_Array1OfShape anArgs (1, aLen);
3164   TCollection_AsciiString anEntry, aPartsDescr, anAsciiList;
3165
3166   for (argi = 1; argi <= aLen; argi++) {
3167     Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast(theParts->Value(argi));
3168     Handle(GEOM_Function) aRef = anObj->GetLastFunction();
3169     if (aRef.IsNull()) return NULL;
3170
3171     TopoDS_Shape anArg = aRef->GetValue();
3172     if (anArg.IsNull()) {
3173       SetErrorCode("Null shape is given as argument");
3174       return NULL;
3175     }
3176     anArgs(argi) = anArg;
3177
3178     // For Python command
3179     TDF_Tool::Entry(anObj->GetEntry(), anEntry);
3180     aPartsDescr += anEntry + ", ";
3181   }
3182
3183   //Get the Blocks
3184   try {
3185 #if OCC_VERSION_LARGE > 0x06010000
3186     OCC_CATCH_SIGNALS;
3187 #endif
3188     TopTools_MapOfShape mapShape;
3189     Standard_Integer nbSolids = 0;
3190     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
3191     for (; exp.More(); exp.Next()) {
3192       if (mapShape.Add(exp.Current())) {
3193         nbSolids++;
3194       }
3195     }
3196
3197     mapShape.Clear();
3198     Standard_Integer ind = 1;
3199     TopTools_Array1OfShape aSolids (1, nbSolids);
3200     TColStd_Array1OfInteger aNbParts (1, nbSolids);
3201     for (exp.Init(aBlockOrComp, TopAbs_SOLID); exp.More(); exp.Next(), ind++) {
3202       if (mapShape.Add(exp.Current())) {
3203         TopoDS_Shape aSolid = exp.Current();
3204         aSolids(ind) = aSolid;
3205         aNbParts(ind) = 0;
3206
3207         // 2. Define quantity of parts, contained in each solid
3208         TopTools_IndexedMapOfShape aSubShapes;
3209         TopExp::MapShapes(aSolid, aSubShapes);
3210         for (argi = 1; argi <= aLen; argi++) {
3211           if (aSubShapes.Contains(anArgs(argi))) {
3212             aNbParts(ind)++;
3213           }
3214         }
3215       }
3216     }
3217
3218     // 3. Define solid, containing maximum quantity of parts
3219     Standard_Integer maxNb = 0, nbFound = 0;
3220     for (ind = 1; ind <= nbSolids; ind++) {
3221       if (aNbParts(ind) > maxNb) {
3222         maxNb = aNbParts(ind);
3223         nbFound = 1;
3224       } else if (aNbParts(ind) == maxNb) {
3225         nbFound++;
3226       } else {
3227       }
3228     }
3229     if (nbFound == 0) {
3230       SetErrorCode("There are no blocks, containing the given parts");
3231       return NULL;
3232     }
3233
3234     // Map shapes
3235     TopTools_IndexedMapOfShape anIndices;
3236     TopExp::MapShapes(aBlockOrComp, anIndices);
3237     Handle(TColStd_HArray1OfInteger) anArray;
3238
3239     for (ind = 1; ind <= nbSolids; ind++) {
3240       if (aNbParts(ind) == maxNb) {
3241         anArray = new TColStd_HArray1OfInteger(1,1);
3242         anArray->SetValue(1, anIndices.FindIndex(aSolids(ind)));
3243         anObj = GetEngine()->AddSubShape(theCompound, anArray);
3244         aBlocks->Append(anObj);
3245
3246         // For Python command
3247         TDF_Tool::Entry(anObj->GetEntry(), anEntry);
3248         anAsciiList += anEntry + ", ";
3249         if (aFunction.IsNull())
3250           aFunction = anObj->GetLastFunction();
3251       }
3252     }
3253   }
3254   catch (Standard_Failure) {
3255     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3256     SetErrorCode(aFail->GetMessageString());
3257     return NULL;
3258   }
3259
3260   //Make a Python command
3261   aPartsDescr.Trunc(aPartsDescr.Length() - 2);
3262   anAsciiList.Trunc(anAsciiList.Length() - 2);
3263
3264   GEOM::TPythonDump(aFunction) << "[" << anAsciiList.ToCString()
3265     << "] = geompy.GetBlocksByParts(" << theCompound
3266       << ", [" << aPartsDescr.ToCString() << "])";
3267
3268   SetErrorCode(OK);
3269   return aBlocks;
3270 }
3271
3272 //=============================================================================
3273 /*!
3274  *  MakeMultiTransformation1D
3275  */
3276 //=============================================================================
3277 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation1D
3278                                              (Handle(GEOM_Object)    theObject,
3279                                               const Standard_Integer theDirFace1,
3280                                               const Standard_Integer theDirFace2,
3281                                               const Standard_Integer theNbTimes)
3282 {
3283   SetErrorCode(KO);
3284
3285   if (theObject.IsNull()) return NULL;
3286
3287   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
3288   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be moved
3289
3290   //Add a new Copy object
3291   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
3292
3293   //Add a translate function
3294   Handle(GEOM_Function) aFunction =
3295     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_MULTI_TRANSFORM_1D);
3296
3297   //Check if the function is set correctly
3298   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
3299
3300   GEOMImpl_IBlockTrsf aTI (aFunction);
3301   aTI.SetOriginal(aLastFunction);
3302   aTI.SetFace1U(theDirFace1);
3303   aTI.SetFace2U(theDirFace2);
3304   aTI.SetNbIterU(theNbTimes);
3305
3306   //Compute the transformation
3307   try {
3308 #if OCC_VERSION_LARGE > 0x06010000
3309     OCC_CATCH_SIGNALS;
3310 #endif
3311     if (!GetSolver()->ComputeFunction(aFunction)) {
3312       SetErrorCode("Block driver failed to make multi-transformation");
3313       return NULL;
3314     }
3315   }
3316   catch (Standard_Failure) {
3317     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3318     SetErrorCode(aFail->GetMessageString());
3319     return NULL;
3320   }
3321
3322   //Make a Python command
3323   GEOM::TPythonDump(aFunction) << aCopy << " = geompy.MakeMultiTransformation1D("
3324     << theObject << ", " << theDirFace1 << ", " << theDirFace2 << ", " << theNbTimes << ")";
3325
3326   SetErrorCode(OK);
3327   return aCopy;
3328 }
3329
3330 //=============================================================================
3331 /*!
3332  *  MakeMultiTransformation2D
3333  */
3334 //=============================================================================
3335 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation2D
3336                                              (Handle(GEOM_Object)    theObject,
3337                                               const Standard_Integer theDirFace1U,
3338                                               const Standard_Integer theDirFace2U,
3339                                               const Standard_Integer theNbTimesU,
3340                                               const Standard_Integer theDirFace1V,
3341                                               const Standard_Integer theDirFace2V,
3342                                               const Standard_Integer theNbTimesV)
3343 {
3344   SetErrorCode(KO);
3345
3346   if (theObject.IsNull()) return NULL;
3347
3348   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
3349   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be moved
3350
3351   //Add a new Copy object
3352   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
3353
3354   //Add a translate function
3355   Handle(GEOM_Function) aFunction =
3356     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_MULTI_TRANSFORM_2D);
3357
3358   //Check if the function is set correctly
3359   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
3360
3361   GEOMImpl_IBlockTrsf aTI (aFunction);
3362   aTI.SetOriginal(aLastFunction);
3363   aTI.SetFace1U(theDirFace1U);
3364   aTI.SetFace2U(theDirFace2U);
3365   aTI.SetNbIterU(theNbTimesU);
3366   aTI.SetFace1V(theDirFace1V);
3367   aTI.SetFace2V(theDirFace2V);
3368   aTI.SetNbIterV(theNbTimesV);
3369
3370   //Compute the transformation
3371   try {
3372 #if OCC_VERSION_LARGE > 0x06010000
3373     OCC_CATCH_SIGNALS;
3374 #endif
3375     if (!GetSolver()->ComputeFunction(aFunction)) {
3376       SetErrorCode("Block driver failed to make multi-transformation");
3377       return NULL;
3378     }
3379   }
3380   catch (Standard_Failure) {
3381     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3382     SetErrorCode(aFail->GetMessageString());
3383     return NULL;
3384   }
3385
3386   //Make a Python command
3387   GEOM::TPythonDump(aFunction) << aCopy << " = geompy.MakeMultiTransformation2D("
3388     << theObject << ", " << theDirFace1U << ", " << theDirFace2U << ", " << theNbTimesU
3389       << ", " << theDirFace1V << ", " << theDirFace2V << ", " << theNbTimesV << ")";
3390
3391   SetErrorCode(OK);
3392   return aCopy;
3393 }
3394
3395 //=============================================================================
3396 /*!
3397  *  Propagate
3398  */
3399 //=============================================================================
3400 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::Propagate
3401                                                  (Handle(GEOM_Object) theShape)
3402 {
3403   SetErrorCode(KO);
3404
3405   if (theShape.IsNull()) return NULL;
3406
3407   TopoDS_Shape aShape = theShape->GetValue();
3408   if (aShape.IsNull()) return NULL;
3409
3410   TopTools_IndexedMapOfShape anIndices;
3411   TopExp::MapShapes(aShape, anIndices);
3412
3413   TopTools_IndexedDataMapOfShapeListOfShape MEW;
3414   GEOMImpl_Block6Explorer::MapShapesAndAncestors
3415     (aShape, TopAbs_EDGE, TopAbs_WIRE, MEW);
3416   Standard_Integer ie, nbEdges = MEW.Extent();
3417
3418   // Result
3419   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
3420
3421   TopTools_MapOfShape mapAcceptedEdges;
3422   TCollection_AsciiString aListRes, anEntry;
3423
3424   // Sort shapes in current chain (Mantis issue 21053)
3425   TopTools_DataMapOfShapeListOfShape aMapChains;
3426   TopTools_ListOfShape aFirstInChains;
3427
3428   for (ie = 1; ie <= nbEdges; ie++) {
3429     TopoDS_Shape curE = MEW.FindKey(ie);
3430
3431     if (mapAcceptedEdges.Contains(curE)) continue;
3432
3433     // Build the chain
3434     TopTools_ListOfShape currentChain;
3435     TopTools_ListOfShape listPrevEdges;
3436
3437     currentChain.Append(curE);
3438     listPrevEdges.Append(curE);
3439     mapAcceptedEdges.Add(curE);
3440
3441     // Collect all edges pass by pass
3442     while (listPrevEdges.Extent() > 0) {
3443       // List of edges, added to chain on this cycle pass
3444       TopTools_ListOfShape listCurEdges;
3445
3446       // Find the next portion of edges
3447       TopTools_ListIteratorOfListOfShape itE (listPrevEdges);
3448       for (; itE.More(); itE.Next()) {
3449         TopoDS_Shape anE = itE.Value();
3450
3451         // Iterate on faces, having edge <anE>
3452         TopTools_ListIteratorOfListOfShape itW (MEW.FindFromKey(anE));
3453         for (; itW.More(); itW.Next()) {
3454           TopoDS_Shape aW = itW.Value();
3455           TopoDS_Shape anOppE;
3456
3457           BRepTools_WireExplorer aWE (TopoDS::Wire(aW));
3458           Standard_Integer nb = 1, found = 0;
3459           TopTools_Array1OfShape anEdges (1,4);
3460           for (; aWE.More(); aWE.Next(), nb++) {
3461             if (nb > 4) {
3462               found = 0;
3463               break;
3464             }
3465             anEdges(nb) = aWE.Current();
3466             if (anEdges(nb).IsSame(anE)) found = nb;
3467           }
3468
3469           if (nb == 5 && found > 0) {
3470             // Quadrangle face found, get an opposite edge
3471             Standard_Integer opp = found + 2;
3472             if (opp > 4) opp -= 4;
3473             anOppE = anEdges(opp);
3474
3475             if (!mapAcceptedEdges.Contains(anOppE)) {
3476               // Add found edge to the chain
3477               currentChain.Append(anOppE);
3478               listCurEdges.Append(anOppE);
3479               mapAcceptedEdges.Add(anOppE);
3480             }
3481           } // if (nb == 5 && found > 0)
3482         } // for (; itF.More(); itF.Next())
3483       } // for (; itE.More(); itE.Next())
3484
3485       listPrevEdges = listCurEdges;
3486     } // while (listPrevEdges.Extent() > 0)
3487
3488     // Sort shapes in current chain (Mantis issue 21053)
3489     GEOMUtils::SortShapes(currentChain, Standard_False);
3490     aFirstInChains.Append(currentChain.First());
3491     aMapChains.Bind(currentChain.First(), currentChain);
3492   }
3493
3494   // Sort chains (Mantis issue 21053)
3495   GEOMUtils::SortShapes(aFirstInChains, Standard_False);
3496
3497   // Store sorted chains in the document
3498   TopTools_ListIteratorOfListOfShape aChainsIt (aFirstInChains);
3499   for (; aChainsIt.More(); aChainsIt.Next()) {
3500     TopoDS_Shape aFirstInChain = aChainsIt.Value();
3501     const TopTools_ListOfShape& currentChain = aMapChains.Find(aFirstInChain);
3502
3503     // Store the chain in the document
3504     Handle(TColStd_HArray1OfInteger) anArray =
3505       new TColStd_HArray1OfInteger (1, currentChain.Extent());
3506
3507     // Fill array of sub-shape indices
3508     TopTools_ListIteratorOfListOfShape itSub (currentChain);
3509     for (int index = 1; itSub.More(); itSub.Next(), ++index) {
3510       int id = anIndices.FindIndex(itSub.Value());
3511       anArray->SetValue(index, id);
3512     }
3513
3514     // Add a new group object
3515     Handle(GEOM_Object) aChain = GetEngine()->AddSubShape(theShape, anArray);
3516
3517     // Set a GROUP type
3518     aChain->SetType(GEOM_GROUP);
3519
3520     // Set a sub-shape type
3521     TDF_Label aFreeLabel = aChain->GetFreeLabel();
3522     TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)TopAbs_EDGE);
3523
3524     // Add the chain to the result
3525     aSeq->Append(aChain);
3526
3527     //Make a Python command
3528     TDF_Tool::Entry(aChain->GetEntry(), anEntry);
3529     aListRes += anEntry + ", ";
3530   }
3531
3532   if (aSeq->IsEmpty()) {
3533     SetErrorCode("There are no quadrangle faces in the shape");
3534     return aSeq;
3535   }
3536
3537   aListRes.Trunc(aListRes.Length() - 2);
3538
3539   // The Propagation doesn't change object so no new function is required.
3540   Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
3541
3542   // Make a Python command
3543   GEOM::TPythonDump(aFunction, /*append=*/true)
3544     << "[" << aListRes.ToCString() << "] = geompy.Propagate(" << theShape << ")";
3545
3546   SetErrorCode(OK);
3547   return aSeq;
3548 }