src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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