src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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