src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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