src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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