src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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