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