src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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