src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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