src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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