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    theTolerance)
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   if (theTolerance < Precision::Confusion()) {
1577     theTolerance == Precision::Confusion();
1578   }
1579
1580   // Compute the result
1581   try {
1582 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
1583     OCC_CATCH_SIGNALS;
1584 #endif
1585     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
1586
1587     TopTools_MapOfShape mapShape;
1588     Standard_Integer nbEdges = 0;
1589     TopExp_Explorer exp (aBlockOrComp, TopAbs_ShapeEnum(theShapeType));
1590     for (; exp.More(); exp.Next()) {
1591       if (mapShape.Add(exp.Current())) {
1592         nbEdges++;
1593       }
1594     }
1595
1596     if (nbEdges == 0) {
1597       SetErrorCode("Given shape contains no subshapes of requested type");
1598       return NULL;
1599     }
1600
1601     // Calculate distances and find min
1602     mapShape.Clear();
1603     Standard_Integer ind = 1;
1604     Standard_Real aMinDist = RealLast();
1605     TopTools_Array1OfShape anEdges (1, nbEdges);
1606     TColStd_Array1OfReal aDistances (1, nbEdges);
1607     for (exp.Init(aBlockOrComp, TopAbs_ShapeEnum(theShapeType)); exp.More(); exp.Next()) {
1608       if (mapShape.Add(exp.Current())) {
1609         TopoDS_Shape anEdge = exp.Current();
1610         anEdges(ind) = anEdge;
1611
1612         BRepExtrema_DistShapeShape aDistTool (aVert, anEdges(ind));
1613         if (!aDistTool.IsDone()) {
1614           SetErrorCode("Can not find a distance from the given point to one of subshapes");
1615           return NULL;
1616         }
1617         aDistances(ind) = aDistTool.Value();
1618         if (aDistances(ind) < aMinDist) {
1619           aMinDist = aDistances(ind);
1620         }
1621         ind++;
1622       }
1623     }
1624
1625     if (aMinDist < RealLast()) {
1626       // Collect subshapes with distance < (aMinDist + theTolerance)
1627       int nbSubShapes = 0;
1628       TopTools_Array1OfShape aNearShapes (1, nbEdges);
1629       for (ind = 1; ind <= nbEdges; ind++) {
1630         if (aDistances(ind) < aMinDist + theTolerance) {
1631           nbSubShapes++;
1632           aNearShapes(nbSubShapes) = anEdges(ind);
1633         }
1634       }
1635
1636       // Add subshape
1637       TopTools_IndexedMapOfShape anIndices;
1638       TopExp::MapShapes(aBlockOrComp, anIndices);
1639       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger (1, nbSubShapes);
1640       for (ind = 1; ind <= nbSubShapes; ind++) {
1641         anArray->SetValue(ind, anIndices.FindIndex(aNearShapes(ind)));
1642       }
1643       aResult = GetEngine()->AddSubShape(theShape, anArray);
1644     }
1645   }
1646   catch (Standard_Failure) {
1647     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1648     SetErrorCode(aFail->GetMessageString());
1649     return NULL;
1650   }
1651
1652   if (aResult.IsNull())
1653     return NULL;
1654
1655   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1656
1657   //Make a Python command
1658   GEOM::TPythonDump(aFunction)
1659     << aResult << " = geompy.GetShapesNearPoint(" << theShape << ", " << thePoint
1660     << ", " << TopAbs_ShapeEnum(theShapeType) << ", " << theTolerance << ")";
1661
1662   SetErrorCode(OK);
1663   return aResult;
1664 }
1665
1666 //=============================================================================
1667 /*!
1668  *  IsCompoundOfBlocks
1669  */
1670 //=============================================================================
1671 Standard_Boolean GEOMImpl_IBlocksOperations::IsCompoundOfBlocks
1672                                                 (Handle(GEOM_Object)    theCompound,
1673                                                  const Standard_Integer theMinNbFaces,
1674                                                  const Standard_Integer theMaxNbFaces,
1675                                                  Standard_Integer&      theNbBlocks)
1676 {
1677   SetErrorCode(KO);
1678   Standard_Boolean isCompOfBlocks = Standard_False;
1679   theNbBlocks = 0;
1680
1681   if (theCompound.IsNull()) return isCompOfBlocks;
1682   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
1683
1684   //Check
1685   isCompOfBlocks = Standard_True;
1686   try {
1687 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
1688     OCC_CATCH_SIGNALS;
1689 #endif
1690     TopTools_MapOfShape mapShape;
1691     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
1692     for (; exp.More(); exp.Next()) {
1693       if (mapShape.Add(exp.Current())) {
1694         TopoDS_Shape aSolid = exp.Current();
1695
1696         TopTools_MapOfShape mapFaces;
1697         TopExp_Explorer expF (aSolid, TopAbs_FACE);
1698         Standard_Integer nbFaces = 0;
1699         for (; expF.More(); expF.Next()) {
1700           if (mapFaces.Add(expF.Current())) {
1701             nbFaces++;
1702             if (nbFaces > theMaxNbFaces) {
1703               isCompOfBlocks = Standard_False;
1704               break;
1705             }
1706           }
1707         }
1708         if (nbFaces < theMinNbFaces || theMaxNbFaces < nbFaces) {
1709           isCompOfBlocks = Standard_False;
1710         } else {
1711           theNbBlocks++;
1712         }
1713       }
1714     }
1715   }
1716   catch (Standard_Failure) {
1717     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1718     SetErrorCode(aFail->GetMessageString());
1719     return isCompOfBlocks;
1720   }
1721
1722   SetErrorCode(OK);
1723   return isCompOfBlocks;
1724 }
1725
1726 //=============================================================================
1727 /*!
1728  *  Set of functions, used by CheckCompoundOfBlocks() method
1729  */
1730 //=============================================================================
1731 void GEOMImpl_IBlocksOperations::AddBlocksFrom (const TopoDS_Shape&   theShape,
1732                                                 TopTools_ListOfShape& BLO,
1733                                                 TopTools_ListOfShape& NOT,
1734                                                 TopTools_ListOfShape& EXT)
1735 {
1736   TopAbs_ShapeEnum aType = theShape.ShapeType();
1737   switch (aType) {
1738   case TopAbs_COMPOUND:
1739   case TopAbs_COMPSOLID:
1740     {
1741       TopoDS_Iterator It (theShape);
1742       for (; It.More(); It.Next()) {
1743         AddBlocksFrom(It.Value(), BLO, NOT, EXT);
1744       }
1745     }
1746     break;
1747   case TopAbs_SOLID:
1748     {
1749       // Check, if there are seam or degenerated edges
1750       BlockFix_CheckTool aTool;
1751       aTool.SetShape(theShape);
1752       aTool.Perform();
1753       if (aTool.NbPossibleBlocks() > 0) {
1754         EXT.Append(theShape);
1755       } else {
1756         // Count faces and edges in each face to recognize blocks
1757         TopTools_MapOfShape mapFaces;
1758         Standard_Integer nbFaces = 0;
1759         Standard_Boolean hasNonQuadr = Standard_False;
1760         TopExp_Explorer expF (theShape, TopAbs_FACE);
1761
1762         for (; expF.More(); expF.Next()) {
1763           if (mapFaces.Add(expF.Current())) {
1764             nbFaces++;
1765             if (nbFaces > 6) break;
1766
1767             // get wire
1768             TopoDS_Shape aF = expF.Current();
1769             TopExp_Explorer wires (aF, TopAbs_WIRE);
1770             if (!wires.More()) {
1771               // no wire in the face
1772               hasNonQuadr = Standard_True;
1773               break;
1774             }
1775             TopoDS_Shape aWire = wires.Current();
1776             wires.Next();
1777             if (wires.More()) {
1778               // multiple wires in the face
1779               hasNonQuadr = Standard_True;
1780               break;
1781             }
1782
1783             // Check number of edges in the face
1784             Standard_Integer nbEdges = 0;
1785             TopTools_MapOfShape mapEdges;
1786             TopExp_Explorer expW (aWire, TopAbs_EDGE);
1787             for (; expW.More(); expW.Next()) {
1788               if (mapEdges.Add(expW.Current())) {
1789                 nbEdges++;
1790                 if (nbEdges > 4) break;
1791               }
1792             }
1793             if (nbEdges != 4) {
1794               hasNonQuadr = Standard_True;
1795             }
1796           }
1797         }
1798
1799         if (nbFaces == 6 && !hasNonQuadr) {
1800           BLO.Append(theShape);
1801         } else {
1802           NOT.Append(theShape);
1803         }
1804       }
1805     }
1806     break;
1807   default:
1808     NOT.Append(theShape);
1809   }
1810 }
1811
1812 void AddBlocksFromOld (const TopoDS_Shape&   theShape,
1813                        TopTools_ListOfShape& BLO,
1814                        TopTools_ListOfShape& NOT,
1815                        TopTools_ListOfShape& DEG,
1816                        TopTools_ListOfShape& SEA)
1817 {
1818   TopAbs_ShapeEnum aType = theShape.ShapeType();
1819   switch (aType) {
1820   case TopAbs_COMPOUND:
1821   case TopAbs_COMPSOLID:
1822     {
1823       TopoDS_Iterator It (theShape);
1824       for (; It.More(); It.Next()) {
1825         AddBlocksFromOld(It.Value(), BLO, NOT, DEG, SEA);
1826       }
1827     }
1828     break;
1829   case TopAbs_SOLID:
1830     {
1831       TopTools_MapOfShape mapFaces;
1832       TopExp_Explorer expF (theShape, TopAbs_FACE);
1833       Standard_Integer nbFaces = 0;
1834       Standard_Boolean hasNonQuadr = Standard_False;
1835       Standard_Boolean hasDegenerated = Standard_False;
1836       Standard_Boolean hasSeam = Standard_False;
1837       for (; expF.More(); expF.Next()) {
1838         if (mapFaces.Add(expF.Current())) {
1839           nbFaces++;
1840           if (nbFaces > 6) break;
1841
1842           // Check number of edges in the face
1843           Standard_Integer nbEdges = 0;
1844           TopTools_MapOfShape mapEdges;
1845
1846           // get wire
1847           TopoDS_Shape aF = expF.Current();
1848           TopExp_Explorer wires (aF, TopAbs_WIRE);
1849           if (!wires.More()) {
1850             // no wire in the face
1851             hasNonQuadr = Standard_True;
1852             break;
1853           }
1854           TopoDS_Shape aWire = wires.Current();
1855           wires.Next();
1856           if (wires.More()) {
1857             // multiple wires in the face
1858             hasNonQuadr = Standard_True;
1859             break;
1860           }
1861
1862           // iterate on wire
1863           BRepTools_WireExplorer aWE (TopoDS::Wire(aWire), TopoDS::Face(aF));
1864           for (; aWE.More(); aWE.Next(), nbEdges++) {
1865             if (BRep_Tool::Degenerated(aWE.Current())) {
1866               // degenerated edge found
1867               hasDegenerated = Standard_True;
1868 //              break;
1869             }
1870             if (mapEdges.Contains(aWE.Current())) {
1871               // seam edge found
1872               hasSeam = Standard_True;
1873 //              break;
1874             }
1875             mapEdges.Add(aWE.Current());
1876           }
1877           if (nbEdges != 4) {
1878             hasNonQuadr = Standard_True;
1879           }
1880         }
1881       }
1882       if (nbFaces == 6) {
1883         if (hasDegenerated || hasSeam) {
1884           if (hasDegenerated) {
1885             DEG.Append(theShape);
1886           }
1887           if (hasSeam) {
1888             SEA.Append(theShape);
1889           }
1890         } else if (hasNonQuadr) {
1891           NOT.Append(theShape);
1892         } else {
1893           BLO.Append(theShape);
1894         }
1895       } else {
1896         NOT.Append(theShape);
1897       }
1898     }
1899     break;
1900   default:
1901     NOT.Append(theShape);
1902   }
1903 }
1904
1905 #define REL_NOT_CONNECTED 0
1906 #define REL_OK            1
1907 #define REL_NOT_GLUED     2
1908 #define REL_COLLISION_VV  3
1909 #define REL_COLLISION_FF  4
1910 #define REL_COLLISION_EE  5
1911 #define REL_UNKNOWN       6
1912
1913 Standard_Integer BlocksRelation (const TopoDS_Shape& theBlock1,
1914                                  const TopoDS_Shape& theBlock2)
1915 {
1916   // Compare bounding boxes before calling BRepExtrema_DistShapeShape
1917   Standard_Real Xmin1, Ymin1, Zmin1, Xmax1, Ymax1, Zmax1;
1918   Standard_Real Xmin2, Ymin2, Zmin2, Xmax2, Ymax2, Zmax2;
1919   Bnd_Box B1, B2;
1920   BRepBndLib::Add(theBlock1, B1);
1921   BRepBndLib::Add(theBlock2, B2);
1922   B1.Get(Xmin1, Ymin1, Zmin1, Xmax1, Ymax1, Zmax1);
1923   B2.Get(Xmin2, Ymin2, Zmin2, Xmax2, Ymax2, Zmax2);
1924   if (Xmax2 < Xmin1 || Xmax1 < Xmin2 ||
1925       Ymax2 < Ymin1 || Ymax1 < Ymin2 ||
1926       Zmax2 < Zmin1 || Zmax1 < Zmin2) {
1927     return REL_NOT_CONNECTED;
1928   }
1929
1930   BRepExtrema_DistShapeShape dst (theBlock1, theBlock2);
1931   if (!dst.IsDone()) {
1932     return REL_UNKNOWN;
1933   }
1934
1935   if (dst.Value() > Precision::Confusion()) {
1936     return REL_NOT_CONNECTED;
1937   }
1938
1939   if (dst.InnerSolution()) {
1940     return REL_COLLISION_VV;
1941   }
1942
1943   Standard_Integer nbSol = dst.NbSolution();
1944   Standard_Integer relation = REL_OK;
1945   Standard_Integer nbVerts = 0;
1946   Standard_Integer nbEdges = 0;
1947   Standard_Integer sol = 1;
1948   for (; sol <= nbSol; sol++) {
1949     BRepExtrema_SupportType supp1 = dst.SupportTypeShape1(sol);
1950     BRepExtrema_SupportType supp2 = dst.SupportTypeShape2(sol);
1951     if (supp1 == BRepExtrema_IsVertex && supp2 == BRepExtrema_IsVertex) {
1952       nbVerts++;
1953     } else if (supp1 == BRepExtrema_IsInFace || supp2 == BRepExtrema_IsInFace) {
1954       return REL_COLLISION_FF;
1955     } else if (supp1 == BRepExtrema_IsOnEdge && supp2 == BRepExtrema_IsOnEdge) {
1956       nbEdges++;
1957     } else if ((supp1 == BRepExtrema_IsOnEdge && supp2 == BRepExtrema_IsVertex) ||
1958                (supp2 == BRepExtrema_IsOnEdge && supp1 == BRepExtrema_IsVertex)) {
1959       relation = REL_COLLISION_EE;
1960     } else {
1961     }
1962   }
1963
1964   if (relation != REL_OK) {
1965     return relation;
1966   }
1967
1968   TColStd_Array1OfInteger vertSol (1, nbVerts);
1969   TopTools_Array1OfShape V1 (1, nbVerts);
1970   TopTools_Array1OfShape V2 (1, nbVerts);
1971   Standard_Integer ivs = 0;
1972   for (sol = 1; sol <= nbSol; sol++) {
1973     if (dst.SupportTypeShape1(sol) == BRepExtrema_IsVertex &&
1974         dst.SupportTypeShape2(sol) == BRepExtrema_IsVertex) {
1975       TopoDS_Vertex Vcur = TopoDS::Vertex(dst.SupportOnShape1(sol));
1976       // Check, that this vertex is far enough from other solution vertices.
1977       Standard_Integer ii = 1;
1978       for (; ii <= ivs; ii++) {
1979         if (BRepTools::Compare(TopoDS::Vertex(V1(ii)), Vcur)) {
1980           continue;
1981         }
1982       }
1983       ivs++;
1984       vertSol(ivs) = sol;
1985       V1(ivs) = Vcur;
1986       V2(ivs) = dst.SupportOnShape2(sol);
1987     }
1988   }
1989
1990   // As we deal only with quadrangles,
1991   // 2, 3 or 4 vertex solutions can be found.
1992   if (ivs <= 1) {
1993     if (nbEdges > 0) {
1994       return REL_COLLISION_FF;
1995     }
1996     return REL_NOT_CONNECTED;
1997   }
1998   if (ivs > 4) {
1999     return REL_UNKNOWN;
2000   }
2001
2002   // Check sharing of coincident entities.
2003   if (ivs == 2 || ivs == 3) {
2004     // Map vertices and edges of the blocks
2005     TopTools_IndexedDataMapOfShapeListOfShape MVE1, MVE2;
2006     GEOMImpl_Block6Explorer::MapShapesAndAncestors
2007       (theBlock1, TopAbs_VERTEX, TopAbs_EDGE, MVE1);
2008     GEOMImpl_Block6Explorer::MapShapesAndAncestors
2009       (theBlock2, TopAbs_VERTEX, TopAbs_EDGE, MVE2);
2010
2011     if (ivs == 2) {
2012       // Find common edge
2013       TopoDS_Shape anEdge1, anEdge2;
2014       GEOMImpl_Block6Explorer::FindEdge(anEdge1, V1(1), V1(2), MVE1);
2015       if (anEdge1.IsNull()) return REL_UNKNOWN;
2016
2017       GEOMImpl_Block6Explorer::FindEdge(anEdge2, V2(1), V2(2), MVE2);
2018       if (anEdge2.IsNull()) return REL_UNKNOWN;
2019
2020       if (!anEdge1.IsSame(anEdge2)) return REL_NOT_GLUED;
2021
2022     } else { // ivs == 3
2023       // Find common edges
2024       Standard_Integer e1_v1 = 1;
2025       Standard_Integer e1_v2 = 2;
2026       Standard_Integer e2_v1 = 3;
2027       Standard_Integer e2_v2 = 1;
2028
2029       TopoDS_Shape anEdge11, anEdge12;
2030       GEOMImpl_Block6Explorer::FindEdge(anEdge11, V1(e1_v1), V1(e1_v2), MVE1);
2031       if (anEdge11.IsNull()) {
2032         e1_v2 = 3;
2033         e2_v1 = 2;
2034         GEOMImpl_Block6Explorer::FindEdge(anEdge11, V1(e1_v1), V1(e1_v2), MVE1);
2035         if (anEdge11.IsNull()) return REL_UNKNOWN;
2036       }
2037       GEOMImpl_Block6Explorer::FindEdge(anEdge12, V1(e2_v1), V1(e2_v2), MVE1);
2038       if (anEdge12.IsNull()) {
2039         e2_v2 = 5 - e2_v1;
2040         GEOMImpl_Block6Explorer::FindEdge(anEdge12, V1(e2_v1), V1(e2_v2), MVE1);
2041         if (anEdge12.IsNull()) return REL_UNKNOWN;
2042       }
2043
2044       TopoDS_Shape anEdge21, anEdge22;
2045       GEOMImpl_Block6Explorer::FindEdge(anEdge21, V2(e1_v1), V2(e1_v2), MVE2);
2046       if (anEdge21.IsNull()) return REL_UNKNOWN;
2047       GEOMImpl_Block6Explorer::FindEdge(anEdge22, V2(e2_v1), V2(e2_v2), MVE2);
2048       if (anEdge22.IsNull()) return REL_UNKNOWN;
2049
2050       // Check of edges coincidence (with some precision) have to be done here
2051       // if (!anEdge11.IsEqual(anEdge21)) return REL_UNKNOWN;
2052       // if (!anEdge12.IsEqual(anEdge22)) return REL_UNKNOWN;
2053
2054       // Check of edges sharing
2055       if (!anEdge11.IsSame(anEdge21)) return REL_NOT_GLUED;
2056       if (!anEdge12.IsSame(anEdge22)) return REL_NOT_GLUED;
2057     }
2058   }
2059
2060   if (ivs == 4) {
2061     // Map vertices and faces of the blocks
2062     TopTools_IndexedDataMapOfShapeListOfShape MVF1, MVF2;
2063     GEOMImpl_Block6Explorer::MapShapesAndAncestors
2064       (theBlock1, TopAbs_VERTEX, TopAbs_FACE, MVF1);
2065     GEOMImpl_Block6Explorer::MapShapesAndAncestors
2066       (theBlock2, TopAbs_VERTEX, TopAbs_FACE, MVF2);
2067
2068     TopoDS_Shape aFace1, aFace2;
2069     GEOMImpl_Block6Explorer::FindFace(aFace1, V1(1), V1(2), V1(3), V1(4), MVF1);
2070     if (aFace1.IsNull()) return REL_UNKNOWN;
2071     GEOMImpl_Block6Explorer::FindFace(aFace2, V2(1), V2(2), V2(3), V2(4), MVF2);
2072     if (aFace2.IsNull()) return REL_UNKNOWN;
2073
2074     // Check of faces coincidence (with some precision) have to be done here
2075     // if (!aFace1.IsEqual(aFace2)) return REL_UNKNOWN;
2076
2077     // Check of faces sharing
2078     if (!aFace1.IsSame(aFace2)) return REL_NOT_GLUED;
2079   }
2080
2081   return REL_OK;
2082 }
2083
2084 void FindConnected (const Standard_Integer         theBlockIndex,
2085                     const TColStd_Array2OfInteger& theRelations,
2086                     TColStd_MapOfInteger&          theProcessedMap,
2087                     TColStd_MapOfInteger&          theConnectedMap)
2088 {
2089   theConnectedMap.Add(theBlockIndex);
2090   theProcessedMap.Add(theBlockIndex);
2091
2092   Standard_Integer nbBlocks = theRelations.ColLength();
2093   Standard_Integer col = 1;
2094   for (; col <= nbBlocks; col++) {
2095     if (theRelations(theBlockIndex, col) == REL_OK ||
2096         theRelations(theBlockIndex, col) == REL_NOT_GLUED) {
2097       if (!theProcessedMap.Contains(col)) {
2098         FindConnected(col, theRelations, theProcessedMap, theConnectedMap);
2099       }
2100     }
2101   }
2102 }
2103
2104 Standard_Boolean HasAnyConnection (const Standard_Integer         theBlockIndex,
2105                                    const TColStd_MapOfInteger&    theWith,
2106                                    const TColStd_Array2OfInteger& theRelations,
2107                                    TColStd_MapOfInteger&          theProcessedMap)
2108 {
2109   theProcessedMap.Add(theBlockIndex);
2110
2111   Standard_Integer nbBlocks = theRelations.ColLength();
2112   Standard_Integer col = 1;
2113   for (; col <= nbBlocks; col++) {
2114     if (theRelations(theBlockIndex, col) != REL_NOT_CONNECTED) {
2115       if (!theProcessedMap.Contains(col)) {
2116         if (theWith.Contains(col))
2117           return Standard_True;
2118         if (HasAnyConnection(col, theWith, theRelations, theProcessedMap))
2119           return Standard_True;
2120       }
2121     }
2122   }
2123
2124   return Standard_False;
2125 }
2126
2127 //=============================================================================
2128 /*!
2129  *  CheckCompoundOfBlocksOld
2130  */
2131 //=============================================================================
2132 Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocksOld
2133                                                 (Handle(GEOM_Object) theCompound,
2134                                                  std::list<BCError>&      theErrors)
2135 {
2136   SetErrorCode(KO);
2137
2138   if (theCompound.IsNull()) return Standard_False;
2139   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2140
2141   Standard_Boolean isCompOfBlocks = Standard_True;
2142
2143   // Map sub-shapes and their indices
2144   TopTools_IndexedMapOfShape anIndices;
2145   TopExp::MapShapes(aBlockOrComp, anIndices);
2146
2147   // 1. Report non-blocks
2148   TopTools_ListOfShape NOT; // Not blocks
2149   TopTools_ListOfShape DEG; // Hexahedral solids, having degenerated edges
2150   TopTools_ListOfShape SEA; // Hexahedral solids, having seam edges
2151   TopTools_ListOfShape BLO; // All blocks from the given compound
2152   AddBlocksFromOld(aBlockOrComp, BLO, NOT, DEG, SEA);
2153
2154   if (NOT.Extent() > 0) {
2155     isCompOfBlocks = Standard_False;
2156     BCError anErr;
2157     anErr.error = NOT_BLOCK;
2158     TopTools_ListIteratorOfListOfShape it (NOT);
2159     for (; it.More(); it.Next()) {
2160       anErr.incriminated.push_back(anIndices.FindIndex(it.Value()));
2161     }
2162     theErrors.push_back(anErr);
2163   }
2164
2165   if (DEG.Extent() > 0 || SEA.Extent() > 0) {
2166     isCompOfBlocks = Standard_False;
2167     BCError anErr;
2168     anErr.error = EXTRA_EDGE;
2169
2170     TopTools_ListIteratorOfListOfShape itDEG (DEG);
2171     for (; itDEG.More(); itDEG.Next()) {
2172       anErr.incriminated.push_back(anIndices.FindIndex(itDEG.Value()));
2173     }
2174
2175     TopTools_ListIteratorOfListOfShape itSEA (SEA);
2176     for (; itSEA.More(); itSEA.Next()) {
2177       anErr.incriminated.push_back(anIndices.FindIndex(itSEA.Value()));
2178     }
2179
2180     theErrors.push_back(anErr);
2181   }
2182
2183   Standard_Integer nbBlocks = BLO.Extent();
2184   if (nbBlocks == 0) {
2185     isCompOfBlocks = Standard_False;
2186     SetErrorCode(OK);
2187     return isCompOfBlocks;
2188   }
2189   if (nbBlocks == 1) {
2190     SetErrorCode(OK);
2191     return isCompOfBlocks;
2192   }
2193
2194   // Convert list of blocks into array for easy and fast access
2195   Standard_Integer ibl = 1;
2196   TopTools_Array1OfShape aBlocks (1, nbBlocks);
2197   TopTools_ListIteratorOfListOfShape BLOit (BLO);
2198   for (; BLOit.More(); BLOit.Next(), ibl++) {
2199     aBlocks.SetValue(ibl, BLOit.Value());
2200   }
2201
2202   // 2. Find relations between all blocks,
2203   //    report connection errors (NOT_GLUED and INVALID_CONNECTION)
2204   TColStd_Array2OfInteger aRelations (1, nbBlocks, 1, nbBlocks);
2205   aRelations.Init(REL_NOT_CONNECTED);
2206
2207   Standard_Integer row = 1;
2208   for (row = 1; row <= nbBlocks; row++) {
2209     TopoDS_Shape aBlock = aBlocks.Value(row);
2210
2211     Standard_Integer col = row + 1;
2212     for (; col <= nbBlocks; col++) {
2213       Standard_Integer aRel = BlocksRelation(aBlock, aBlocks.Value(col));
2214       if (aRel != REL_NOT_CONNECTED) {
2215         aRelations.SetValue(row, col, aRel);
2216         aRelations.SetValue(col, row, aRel);
2217         if (aRel == REL_NOT_GLUED) {
2218           // report connection error
2219           isCompOfBlocks = Standard_False;
2220           BCError anErr;
2221           anErr.error = NOT_GLUED;
2222           anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(row)));
2223           anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(col)));
2224           theErrors.push_back(anErr);
2225         } else if (aRel == REL_COLLISION_VV ||
2226                    aRel == REL_COLLISION_FF ||
2227                    aRel == REL_COLLISION_EE ||
2228                    aRel == REL_UNKNOWN) {
2229           // report connection error
2230           isCompOfBlocks = Standard_False;
2231           BCError anErr;
2232           anErr.error = INVALID_CONNECTION;
2233           anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(row)));
2234           anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(col)));
2235           theErrors.push_back(anErr);
2236         } else {
2237         }
2238       }
2239     }
2240   }
2241
2242   // 3. Find largest set of connected (good connection or not glued) blocks
2243   TColStd_MapOfInteger aProcessedMap;
2244   TColStd_MapOfInteger aLargestSet;
2245   TColStd_MapOfInteger aCurrentSet;
2246   for (ibl = 1; ibl <= nbBlocks; ibl++) {
2247     if (!aProcessedMap.Contains(ibl)) {
2248       aCurrentSet.Clear();
2249       FindConnected(ibl, aRelations, aProcessedMap, aCurrentSet);
2250       if (aCurrentSet.Extent() > aLargestSet.Extent()) {
2251         aLargestSet = aCurrentSet;
2252       }
2253     }
2254   }
2255
2256   // 4. Report all blocks, isolated from <aLargestSet>
2257   BCError anErr;
2258   anErr.error = NOT_CONNECTED;
2259   Standard_Boolean hasIsolated = Standard_False;
2260   for (ibl = 1; ibl <= nbBlocks; ibl++) {
2261     if (!aLargestSet.Contains(ibl)) {
2262       aProcessedMap.Clear();
2263       if (!HasAnyConnection(ibl, aLargestSet, aRelations, aProcessedMap)) {
2264         // report connection absence
2265         hasIsolated = Standard_True;
2266         anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(ibl)));
2267       }
2268     }
2269   }
2270   if (hasIsolated) {
2271     isCompOfBlocks = Standard_False;
2272     theErrors.push_back(anErr);
2273   }
2274
2275   SetErrorCode(OK);
2276   return isCompOfBlocks;
2277 }
2278
2279 //=============================================================================
2280 /*!
2281  *  PrintBCErrors
2282  */
2283 //=============================================================================
2284 TCollection_AsciiString GEOMImpl_IBlocksOperations::PrintBCErrors
2285                                                 (Handle(GEOM_Object)  theCompound,
2286                                                  const std::list<BCError>& theErrors)
2287 {
2288   TCollection_AsciiString aDescr;
2289
2290   std::list<BCError>::const_iterator errIt = theErrors.begin();
2291   int i = 0;
2292   for (; errIt != theErrors.end(); i++, errIt++) {
2293     BCError errStruct = *errIt;
2294
2295     switch (errStruct.error) {
2296     case NOT_BLOCK:
2297       aDescr += "\n\tNot a Blocks: ";
2298       break;
2299     case EXTRA_EDGE:
2300       aDescr += "\n\tHexahedral solids with degenerated and/or seam edges: ";
2301       break;
2302     case INVALID_CONNECTION:
2303       aDescr += "\n\tInvalid connection between two blocks: ";
2304       break;
2305     case NOT_CONNECTED:
2306       aDescr += "\n\tBlocks, not connected with main body: ";
2307       break;
2308     case NOT_GLUED:
2309       aDescr += "\n\tNot glued blocks: ";
2310       break;
2311     default:
2312       break;
2313     }
2314
2315     std::list<int> sshList = errStruct.incriminated;
2316     std::list<int>::iterator sshIt = sshList.begin();
2317     int jj = 0;
2318     for (; sshIt != sshList.end(); jj++, sshIt++) {
2319       if (jj > 0)
2320         aDescr += ", ";
2321       aDescr += TCollection_AsciiString(*sshIt);
2322     }
2323   }
2324
2325   return aDescr;
2326 }
2327
2328 //=============================================================================
2329 /*!
2330  *  CheckCompoundOfBlocks
2331  */
2332 //=============================================================================
2333 Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
2334                                               (Handle(GEOM_Object) theCompound,
2335                                                std::list<BCError>& theErrors)
2336 {
2337   SetErrorCode(KO);
2338
2339   if (theCompound.IsNull()) return Standard_False;
2340   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2341
2342   Standard_Boolean isCompOfBlocks = Standard_True;
2343
2344   // Map sub-shapes and their indices
2345   TopTools_IndexedMapOfShape anIndices;
2346   TopExp::MapShapes(aBlockOrComp, anIndices);
2347
2348   // 1. Separate blocks from non-blocks
2349   TopTools_ListOfShape NOT; // Not blocks
2350   TopTools_ListOfShape EXT; // Hexahedral solids, having degenerated and/or seam edges
2351   TopTools_ListOfShape BLO; // All blocks from the given compound
2352   AddBlocksFrom(aBlockOrComp, BLO, NOT, EXT);
2353
2354   // Report non-blocks
2355   if (NOT.Extent() > 0) {
2356     isCompOfBlocks = Standard_False;
2357     BCError anErr;
2358     anErr.error = NOT_BLOCK;
2359     TopTools_ListIteratorOfListOfShape it (NOT);
2360     for (; it.More(); it.Next()) {
2361       anErr.incriminated.push_back(anIndices.FindIndex(it.Value()));
2362     }
2363     theErrors.push_back(anErr);
2364   }
2365
2366   // Report solids, having degenerated and/or seam edges
2367   if (EXT.Extent() > 0) {
2368     isCompOfBlocks = Standard_False;
2369     BCError anErr;
2370     anErr.error = EXTRA_EDGE;
2371     TopTools_ListIteratorOfListOfShape it (EXT);
2372     for (; it.More(); it.Next()) {
2373       anErr.incriminated.push_back(anIndices.FindIndex(it.Value()));
2374     }
2375     theErrors.push_back(anErr);
2376   }
2377
2378   Standard_Integer nbBlocks = BLO.Extent();
2379   if (nbBlocks == 0) {
2380     isCompOfBlocks = Standard_False;
2381     SetErrorCode(OK);
2382     return isCompOfBlocks;
2383   }
2384   if (nbBlocks == 1) {
2385     SetErrorCode(OK);
2386     return isCompOfBlocks;
2387   }
2388
2389   // Prepare data for 2. and 3.
2390   TColStd_Array2OfInteger aRelations (1, nbBlocks, 1, nbBlocks);
2391   aRelations.Init(REL_NOT_CONNECTED);
2392
2393   TopTools_IndexedMapOfShape mapBlocks;
2394
2395   BRep_Builder BB;
2396   TopoDS_Compound aComp;
2397   BB.MakeCompound(aComp);
2398
2399   TopTools_ListIteratorOfListOfShape BLOit (BLO);
2400   for (; BLOit.More(); BLOit.Next()) {
2401     mapBlocks.Add(BLOit.Value());
2402     BB.Add(aComp, BLOit.Value());
2403   }
2404
2405   // 2. Find glued blocks (having shared faces)
2406   TopTools_IndexedDataMapOfShapeListOfShape mapFaceBlocks;
2407   GEOMImpl_Block6Explorer::MapShapesAndAncestors
2408     (aComp, TopAbs_FACE, TopAbs_SOLID, mapFaceBlocks);
2409
2410   Standard_Integer prevInd = 0, curInd = 0;
2411   Standard_Integer ind = 1, nbFaces = mapFaceBlocks.Extent();
2412   for (; ind <= nbFaces; ind++) {
2413     const TopTools_ListOfShape& aGluedBlocks = mapFaceBlocks.FindFromIndex(ind);
2414     if (aGluedBlocks.Extent() > 1) { // Shared face found
2415       TopTools_ListIteratorOfListOfShape aGluedBlocksIt (aGluedBlocks);
2416       TopoDS_Shape prevBlock, curBlock;
2417       for (; aGluedBlocksIt.More(); aGluedBlocksIt.Next()) {
2418         curBlock = aGluedBlocksIt.Value();
2419         if (!prevBlock.IsNull()) {
2420           prevInd = mapBlocks.FindIndex(prevBlock);
2421           curInd  = mapBlocks.FindIndex(curBlock);
2422           aRelations.SetValue(prevInd, curInd, REL_OK);
2423           aRelations.SetValue(curInd, prevInd, REL_OK);
2424         }
2425         prevBlock = curBlock;
2426       }
2427     }
2428   }
2429
2430   // 3. Find not glued blocks
2431   GEOMAlgo_GlueAnalyser aGD;
2432
2433   aGD.SetShape(aComp);
2434   aGD.SetTolerance(Precision::Confusion());
2435   aGD.SetCheckGeometry(Standard_True);
2436   aGD.Perform();
2437
2438   Standard_Integer iErr, iWrn;
2439   iErr = aGD.ErrorStatus();
2440   if (iErr) {
2441     SetErrorCode("Error in GEOMAlgo_GlueAnalyser");
2442     return isCompOfBlocks;
2443   }
2444   iWrn = aGD.WarningStatus();
2445   if (iWrn) {
2446     MESSAGE("Warning in GEOMAlgo_GlueAnalyser");
2447   }
2448
2449   // Report not glued blocks
2450   if (aGD.HasSolidsToGlue()) {
2451     isCompOfBlocks = Standard_False;
2452     Standard_Integer aSx1Ind, aSx2Ind;
2453
2454     const GEOMAlgo_ListOfCoupleOfShapes& aLCS = aGD.SolidsToGlue();
2455     GEOMAlgo_ListIteratorOfListOfCoupleOfShapes aItCS (aLCS);
2456     for (; aItCS.More(); aItCS.Next()) {
2457       const GEOMAlgo_CoupleOfShapes& aCS = aItCS.Value();
2458       const TopoDS_Shape& aSx1 = aCS.Shape1();
2459       const TopoDS_Shape& aSx2 = aCS.Shape2();
2460
2461       aSx1Ind = mapBlocks.FindIndex(aSx1);
2462       aSx2Ind = mapBlocks.FindIndex(aSx2);
2463       aRelations.SetValue(aSx1Ind, aSx2Ind, NOT_GLUED);
2464       aRelations.SetValue(aSx2Ind, aSx1Ind, NOT_GLUED);
2465
2466       BCError anErr;
2467       anErr.error = NOT_GLUED;
2468       anErr.incriminated.push_back(anIndices.FindIndex(aSx1));
2469       anErr.incriminated.push_back(anIndices.FindIndex(aSx2));
2470       theErrors.push_back(anErr);
2471     }
2472   }
2473
2474   // 4. Find largest set of connected (good connection or not glued) blocks
2475   Standard_Integer ibl = 1;
2476   TColStd_MapOfInteger aProcessedMap;
2477   TColStd_MapOfInteger aLargestSet;
2478   TColStd_MapOfInteger aCurrentSet;
2479   for (ibl = 1; ibl <= nbBlocks; ibl++) {
2480     if (!aProcessedMap.Contains(ibl)) {
2481       aCurrentSet.Clear();
2482       FindConnected(ibl, aRelations, aProcessedMap, aCurrentSet);
2483       if (aCurrentSet.Extent() > aLargestSet.Extent()) {
2484         aLargestSet = aCurrentSet;
2485       }
2486     }
2487   }
2488
2489   // 5. Report all blocks, isolated from <aLargestSet>
2490   BCError anErr;
2491   anErr.error = NOT_CONNECTED;
2492   Standard_Boolean hasIsolated = Standard_False;
2493   for (ibl = 1; ibl <= nbBlocks; ibl++) {
2494     if (!aLargestSet.Contains(ibl)) {
2495       aProcessedMap.Clear();
2496       if (!HasAnyConnection(ibl, aLargestSet, aRelations, aProcessedMap)) {
2497         // report connection absence
2498         hasIsolated = Standard_True;
2499         anErr.incriminated.push_back(anIndices.FindIndex(mapBlocks.FindKey(ibl)));
2500       }
2501     }
2502   }
2503   if (hasIsolated) {
2504     isCompOfBlocks = Standard_False;
2505     theErrors.push_back(anErr);
2506   }
2507
2508   SetErrorCode(OK);
2509   return isCompOfBlocks;
2510 }
2511
2512 //=============================================================================
2513 /*!
2514  *  RemoveExtraEdges
2515  */
2516 //=============================================================================
2517 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::RemoveExtraEdges
2518                                      (Handle(GEOM_Object) theObject,
2519                                       const Standard_Integer theOptimumNbFaces)
2520 {
2521   SetErrorCode(KO);
2522
2523   if (theObject.IsNull()) return NULL;
2524
2525   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
2526   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be fixed
2527
2528   //Add a new Copy object
2529   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
2530
2531   //Add a function
2532   Handle(GEOM_Function) aFunction =
2533     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_REMOVE_EXTRA);
2534
2535   //Check if the function is set correctly
2536   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
2537
2538   GEOMImpl_IBlockTrsf aTI (aFunction);
2539   aTI.SetOriginal(aLastFunction);
2540   aTI.SetOptimumNbFaces(theOptimumNbFaces);
2541
2542   //Compute the fixed shape
2543   try {
2544 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
2545     OCC_CATCH_SIGNALS;
2546 #endif
2547     if (!GetSolver()->ComputeFunction(aFunction)) {
2548       SetErrorCode("Block driver failed to remove extra edges of the given shape");
2549       return NULL;
2550     }
2551   }
2552   catch (Standard_Failure) {
2553     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2554     SetErrorCode(aFail->GetMessageString());
2555     return NULL;
2556   }
2557
2558   //Make a Python command
2559   std::string doUnionFaces = (theOptimumNbFaces < 0) ? "False" : "True";
2560   GEOM::TPythonDump(aFunction) << aCopy << " = geompy.RemoveExtraEdges("
2561                                << theObject << ", " << doUnionFaces.data() << ")";
2562
2563   SetErrorCode(OK);
2564   return aCopy;
2565 }
2566
2567 //=============================================================================
2568 /*!
2569  *  CheckAndImprove
2570  */
2571 //=============================================================================
2572 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::CheckAndImprove
2573                                              (Handle(GEOM_Object) theObject)
2574 {
2575   SetErrorCode(KO);
2576
2577   if (theObject.IsNull()) return NULL;
2578
2579   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
2580   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be fixed
2581
2582   //Add a new Copy object
2583   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
2584
2585   //Add a function
2586   Handle(GEOM_Function) aFunction =
2587     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_COMPOUND_IMPROVE);
2588
2589   //Check if the function is set correctly
2590   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
2591
2592   GEOMImpl_IBlockTrsf aTI (aFunction);
2593   aTI.SetOriginal(aLastFunction);
2594
2595   // -1 means do not unite faces on common surface (?except case of seam edge between them?)
2596   //aTI.SetOptimumNbFaces(-1);
2597   aTI.SetOptimumNbFaces(6);
2598
2599   //Compute the fixed shape
2600   try {
2601 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
2602     OCC_CATCH_SIGNALS;
2603 #endif
2604     if (!GetSolver()->ComputeFunction(aFunction)) {
2605       SetErrorCode("Block driver failed to improve the given blocks compound");
2606       return NULL;
2607     }
2608   }
2609   catch (Standard_Failure) {
2610     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2611     SetErrorCode(aFail->GetMessageString());
2612     return NULL;
2613   }
2614
2615   //Make a Python command
2616   GEOM::TPythonDump(aFunction) << aCopy
2617     << " = geompy.CheckAndImprove(" << theObject << ")";
2618
2619   SetErrorCode(OK);
2620   return aCopy;
2621 }
2622
2623 //=============================================================================
2624 /*!
2625  *  ExplodeCompoundOfBlocks
2626  */
2627 //=============================================================================
2628 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::ExplodeCompoundOfBlocks
2629                                                 (Handle(GEOM_Object)    theCompound,
2630                                                  const Standard_Integer theMinNbFaces,
2631                                                  const Standard_Integer theMaxNbFaces)
2632 {
2633   SetErrorCode(KO);
2634
2635   if (theCompound.IsNull()) return NULL;
2636   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2637   if (aBlockOrComp.IsNull()) return NULL;
2638
2639   Handle(TColStd_HSequenceOfTransient) aBlocks = new TColStd_HSequenceOfTransient;
2640   Handle(GEOM_Object) anObj;
2641   Handle(GEOM_Function) aFunction;
2642
2643   TopTools_MapOfShape mapShape;
2644   TCollection_AsciiString anAsciiList, anEntry;
2645
2646   // Map shapes
2647   TopTools_IndexedMapOfShape anIndices;
2648   TopExp::MapShapes(aBlockOrComp, anIndices);
2649   Handle(TColStd_HArray1OfInteger) anArray;
2650
2651   // Explode
2652   try {
2653 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
2654     OCC_CATCH_SIGNALS;
2655 #endif
2656     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
2657     for (; exp.More(); exp.Next()) {
2658       if (mapShape.Add(exp.Current())) {
2659         TopoDS_Shape aSolid = exp.Current();
2660
2661         TopTools_MapOfShape mapFaces;
2662         TopExp_Explorer expF (aSolid, TopAbs_FACE);
2663         Standard_Integer nbFaces = 0;
2664         for (; expF.More(); expF.Next()) {
2665           if (mapFaces.Add(expF.Current())) {
2666             nbFaces++;
2667           }
2668         }
2669
2670         if (theMinNbFaces <= nbFaces && nbFaces <= theMaxNbFaces) {
2671           anArray = new TColStd_HArray1OfInteger(1,1);
2672           anArray->SetValue(1, anIndices.FindIndex(aSolid));
2673           anObj = GetEngine()->AddSubShape(theCompound, anArray);
2674           aBlocks->Append(anObj);
2675
2676           //Make a Python command
2677           TDF_Tool::Entry(anObj->GetEntry(), anEntry);
2678           anAsciiList += anEntry + ", ";
2679         }
2680       }
2681     }
2682   }
2683   catch (Standard_Failure) {
2684     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2685     SetErrorCode(aFail->GetMessageString());
2686     return aBlocks;
2687   }
2688
2689   if (aBlocks->IsEmpty()) {
2690     SetErrorCode("There are no specified blocks in the given shape");
2691     return aBlocks;
2692   }
2693
2694   anAsciiList.Trunc(anAsciiList.Length() - 2);
2695
2696   //The explode doesn't change object so no new function is required.
2697   aFunction = theCompound->GetLastFunction();
2698
2699   //Make a Python command
2700   GEOM::TPythonDump(aFunction, /*append=*/true)
2701     << "[" << anAsciiList.ToCString() << "] = geompy.MakeBlockExplode("
2702     << theCompound << ", " << theMinNbFaces << ", " << theMaxNbFaces << ")";
2703
2704   SetErrorCode(OK);
2705   return aBlocks;
2706 }
2707
2708 //=============================================================================
2709 /*!
2710  *  GetBlockNearPoint
2711  */
2712 //=============================================================================
2713 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockNearPoint
2714                                                 (Handle(GEOM_Object) theCompound,
2715                                                  Handle(GEOM_Object) thePoint)
2716 {
2717   SetErrorCode(KO);
2718
2719   //New object
2720   Handle(GEOM_Object) aResult;
2721
2722   // Arguments
2723   if (theCompound.IsNull() || thePoint.IsNull()) return NULL;
2724
2725   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2726   if (aBlockOrComp.IsNull()) {
2727     SetErrorCode("Compound is null");
2728     return NULL;
2729   }
2730   if (aBlockOrComp.ShapeType() != TopAbs_COMPOUND &&
2731       aBlockOrComp.ShapeType() != TopAbs_COMPSOLID) {
2732     SetErrorCode("Shape to find block in is not a compound");
2733     return NULL;
2734   }
2735
2736   TopoDS_Shape anArg = thePoint->GetValue();
2737   if (anArg.IsNull()) {
2738     SetErrorCode("Point is null");
2739     return NULL;
2740   }
2741   if (anArg.ShapeType() != TopAbs_VERTEX) {
2742     SetErrorCode("Shape for block identification is not a vertex");
2743     return NULL;
2744   }
2745
2746   //Compute the Block value
2747   try {
2748 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
2749     OCC_CATCH_SIGNALS;
2750 #endif
2751     TopoDS_Shape aShape;
2752
2753     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
2754     gp_Pnt aPnt = BRep_Tool::Pnt(aVert);
2755     Standard_Real PX, PY, PZ;
2756     aPnt.Coord(PX, PY, PZ);
2757
2758     // 1. Classify the point relatively each block
2759     Standard_Integer nearest = 2, nbFound = 0;
2760     TopTools_DataMapOfShapeInteger mapShapeDist;
2761     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
2762     for (; exp.More(); exp.Next()) {
2763       TopoDS_Shape aSolid = exp.Current();
2764
2765       if (!mapShapeDist.IsBound(aSolid)) {
2766         Standard_Integer aDistance = 2;
2767
2768         // 1.a. Classify relatively Bounding box
2769         Standard_Real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
2770         Bnd_Box BB;
2771         BRepBndLib::Add(aSolid, BB);
2772         BB.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
2773         if (PX < Xmin || Xmax < PX ||
2774             PY < Ymin || Ymax < PY ||
2775             PZ < Zmin || Zmax < PZ) {
2776           // OUT of bounding box
2777           aDistance = 1;
2778         } else {
2779           // 1.b. Classify relatively the solid itself
2780           BRepClass3d_SolidClassifier SC (aSolid, aPnt, Precision::Confusion());
2781           if (SC.State() == TopAbs_IN) {
2782             aDistance = -1;
2783           } else if (SC.State() == TopAbs_ON) {
2784             aDistance = 0;
2785           } else { // OUT
2786             aDistance = 1;
2787           }
2788         }
2789
2790         if (aDistance < nearest) {
2791           nearest = aDistance;
2792           aShape = aSolid;
2793           nbFound = 1;
2794
2795           // A first found block, containing the point inside, will be returned.
2796           // It is the solution, if there are no intersecting blocks in the compound.
2797           if (nearest == -1) break;
2798
2799         } else if (aDistance == nearest) {
2800           nbFound++;
2801         } else {
2802         }
2803
2804         mapShapeDist.Bind(aSolid, aDistance);
2805       } // if (!mapShapeDist.IsBound(aSolid))
2806     }
2807
2808     // 2. Define block, containing the point or having minimum distance to it
2809     if (nbFound > 1) {
2810       if (nearest == 0) {
2811         // The point is on boundary of some blocks and there are
2812         // no blocks, having the point inside their volume
2813         SetErrorCode("Multiple blocks near the given point are found");
2814         return NULL;
2815
2816       } else if (nearest == 1) {
2817         // The point is outside some blocks and there are
2818         // no blocks, having the point inside or on boundary.
2819         // We will get a nearest block
2820         Standard_Real minDist = RealLast();
2821         TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
2822         for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
2823           if (mapShapeDistIter.Value() == 1) {
2824             TopoDS_Shape aSolid = mapShapeDistIter.Key();
2825             BRepExtrema_DistShapeShape aDistTool (aVert, aSolid);
2826             if (!aDistTool.IsDone()) {
2827               SetErrorCode("Can not find a distance from the given point to one of blocks");
2828               return NULL;
2829             }
2830             Standard_Real aDist = aDistTool.Value();
2831             if (aDist < minDist) {
2832               minDist = aDist;
2833               aShape = aSolid;
2834             }
2835           }
2836         }
2837       } else { // nearest == -1
2838 //        // The point is inside some blocks.
2839 //        // We will get a block with nearest center
2840 //        Standard_Real minDist = RealLast();
2841 //        TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
2842 //        for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
2843 //          if (mapShapeDistIter.Value() == -1) {
2844 //            TopoDS_Shape aSolid = mapShapeDistIter.Key();
2845 //            GProp_GProps aSystem;
2846 //            BRepGProp::VolumeProperties(aSolid, aSystem);
2847 //            gp_Pnt aCenterMass = aSystem.CentreOfMass();
2848 //
2849 //            Standard_Real aDist = aCenterMass.Distance(aPnt);
2850 //            if (aDist < minDist) {
2851 //              minDist = aDist;
2852 //              aShape = aSolid;
2853 //            }
2854 //          }
2855 //        }
2856       }
2857     } // if (nbFound > 1)
2858
2859     if (nbFound == 0) {
2860       SetErrorCode("There are no blocks near the given point");
2861       return NULL;
2862     } else {
2863       TopTools_IndexedMapOfShape anIndices;
2864       TopExp::MapShapes(aBlockOrComp, anIndices);
2865       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
2866       anArray->SetValue(1, anIndices.FindIndex(aShape));
2867       aResult = GetEngine()->AddSubShape(theCompound, anArray);
2868     }
2869   }
2870   catch (Standard_Failure) {
2871     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2872     SetErrorCode(aFail->GetMessageString());
2873     return NULL;
2874   }
2875
2876   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
2877
2878   //Make a Python command
2879   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetBlockNearPoint("
2880                                << theCompound << ", " << thePoint << ")";
2881
2882   SetErrorCode(OK);
2883   return aResult;
2884 }
2885
2886 //=============================================================================
2887 /*!
2888  *  GetBlockByParts
2889  */
2890 //=============================================================================
2891 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockByParts
2892                       (Handle(GEOM_Object)                         theCompound,
2893                        const Handle(TColStd_HSequenceOfTransient)& theParts)
2894 {
2895   SetErrorCode(KO);
2896
2897   Handle(GEOM_Object) aResult;
2898
2899   if (theCompound.IsNull() || theParts.IsNull()) return NULL;
2900   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2901   if (aBlockOrComp.IsNull()) return NULL;
2902
2903   //Get the parts
2904   Standard_Integer argi, aLen = theParts->Length();
2905   TopTools_Array1OfShape anArgs (1, aLen);
2906   TCollection_AsciiString anEntry, aPartsDescr;
2907   for (argi = 1; argi <= aLen; argi++) {
2908     Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast(theParts->Value(argi));
2909     Handle(GEOM_Function) aRef = anObj->GetLastFunction();
2910     if (aRef.IsNull()) return NULL;
2911
2912     TopoDS_Shape anArg = aRef->GetValue();
2913     if (anArg.IsNull()) {
2914       SetErrorCode("Null shape is given as argument");
2915       return NULL;
2916     }
2917     anArgs(argi) = anArg;
2918
2919     // For Python command
2920     TDF_Tool::Entry(anObj->GetEntry(), anEntry);
2921     if (argi > 1) aPartsDescr += ", ";
2922     aPartsDescr += anEntry;
2923   }
2924
2925   //Compute the Block value
2926   try {
2927 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
2928     OCC_CATCH_SIGNALS;
2929 #endif
2930     // 1. Explode compound on solids
2931     TopTools_MapOfShape mapShape;
2932     Standard_Integer nbSolids = 0;
2933     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
2934     for (; exp.More(); exp.Next()) {
2935       if (mapShape.Add(exp.Current())) {
2936         nbSolids++;
2937       }
2938     }
2939
2940     mapShape.Clear();
2941     Standard_Integer ind = 1;
2942     TopTools_Array1OfShape aSolids (1, nbSolids);
2943     TColStd_Array1OfInteger aNbParts (1, nbSolids);
2944     for (exp.Init(aBlockOrComp, TopAbs_SOLID); exp.More(); exp.Next(), ind++) {
2945       if (mapShape.Add(exp.Current())) {
2946         TopoDS_Shape aSolid = exp.Current();
2947         aSolids(ind) = aSolid;
2948         aNbParts(ind) = 0;
2949
2950         // 2. Define quantity of parts, contained in each solid
2951         TopTools_IndexedMapOfShape aSubShapes;
2952         TopExp::MapShapes(aSolid, aSubShapes);
2953         for (argi = 1; argi <= aLen; argi++) {
2954           if (aSubShapes.Contains(anArgs(argi))) {
2955             aNbParts(ind)++;
2956           }
2957         }
2958       }
2959     }
2960
2961     // 3. Define solid, containing maximum quantity of parts
2962     Standard_Integer maxNb = 0, nbFound = 0;
2963     TopoDS_Shape aShape;
2964     for (ind = 1; ind <= nbSolids; ind++) {
2965       if (aNbParts(ind) > maxNb) {
2966         maxNb = aNbParts(ind);
2967         aShape = aSolids(ind);
2968         nbFound = 1;
2969       } else if (aNbParts(ind) == maxNb) {
2970         nbFound++;
2971       } else {
2972       }
2973     }
2974     if (nbFound > 1) {
2975       SetErrorCode("Multiple blocks, containing maximum quantity of the given parts, are found");
2976       return NULL;
2977     } else if (nbFound == 0) {
2978       SetErrorCode("There are no blocks, containing the given parts");
2979       return NULL;
2980     } else {
2981       TopTools_IndexedMapOfShape anIndices;
2982       TopExp::MapShapes(aBlockOrComp, anIndices);
2983       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
2984       anArray->SetValue(1, anIndices.FindIndex(aShape));
2985       aResult = GetEngine()->AddSubShape(theCompound, anArray);
2986     }
2987   } catch (Standard_Failure) {
2988     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2989     SetErrorCode(aFail->GetMessageString());
2990     return NULL;
2991   }
2992
2993   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
2994
2995   //Make a Python command
2996   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetBlockByParts("
2997     << theCompound << ", [" << aPartsDescr.ToCString() << "])";
2998
2999   SetErrorCode(OK);
3000   return aResult;
3001 }
3002
3003 //=============================================================================
3004 /*!
3005  *  GetBlocksByParts
3006  */
3007 //=============================================================================
3008 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::GetBlocksByParts
3009                       (Handle(GEOM_Object)                         theCompound,
3010                        const Handle(TColStd_HSequenceOfTransient)& theParts)
3011 {
3012   SetErrorCode(KO);
3013
3014   if (theCompound.IsNull() || theParts.IsNull()) return NULL;
3015   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
3016   if (aBlockOrComp.IsNull()) return NULL;
3017
3018   Handle(TColStd_HSequenceOfTransient) aBlocks = new TColStd_HSequenceOfTransient;
3019   Handle(GEOM_Object) anObj;
3020   Handle(GEOM_Function) aFunction;
3021
3022   //Get the parts
3023   Standard_Integer argi, aLen = theParts->Length();
3024   TopTools_Array1OfShape anArgs (1, aLen);
3025   TCollection_AsciiString anEntry, aPartsDescr, anAsciiList;
3026
3027   for (argi = 1; argi <= aLen; argi++) {
3028     Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast(theParts->Value(argi));
3029     Handle(GEOM_Function) aRef = anObj->GetLastFunction();
3030     if (aRef.IsNull()) return NULL;
3031
3032     TopoDS_Shape anArg = aRef->GetValue();
3033     if (anArg.IsNull()) {
3034       SetErrorCode("Null shape is given as argument");
3035       return NULL;
3036     }
3037     anArgs(argi) = anArg;
3038
3039     // For Python command
3040     TDF_Tool::Entry(anObj->GetEntry(), anEntry);
3041     aPartsDescr += anEntry + ", ";
3042   }
3043
3044   //Get the Blocks
3045   try {
3046 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
3047     OCC_CATCH_SIGNALS;
3048 #endif
3049     TopTools_MapOfShape mapShape;
3050     Standard_Integer nbSolids = 0;
3051     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
3052     for (; exp.More(); exp.Next()) {
3053       if (mapShape.Add(exp.Current())) {
3054         nbSolids++;
3055       }
3056     }
3057
3058     mapShape.Clear();
3059     Standard_Integer ind = 1;
3060     TopTools_Array1OfShape aSolids (1, nbSolids);
3061     TColStd_Array1OfInteger aNbParts (1, nbSolids);
3062     for (exp.Init(aBlockOrComp, TopAbs_SOLID); exp.More(); exp.Next(), ind++) {
3063       if (mapShape.Add(exp.Current())) {
3064         TopoDS_Shape aSolid = exp.Current();
3065         aSolids(ind) = aSolid;
3066         aNbParts(ind) = 0;
3067
3068         // 2. Define quantity of parts, contained in each solid
3069         TopTools_IndexedMapOfShape aSubShapes;
3070         TopExp::MapShapes(aSolid, aSubShapes);
3071         for (argi = 1; argi <= aLen; argi++) {
3072           if (aSubShapes.Contains(anArgs(argi))) {
3073             aNbParts(ind)++;
3074           }
3075         }
3076       }
3077     }
3078
3079     // 3. Define solid, containing maximum quantity of parts
3080     Standard_Integer maxNb = 0, nbFound = 0;
3081     for (ind = 1; ind <= nbSolids; ind++) {
3082       if (aNbParts(ind) > maxNb) {
3083         maxNb = aNbParts(ind);
3084         nbFound = 1;
3085       } else if (aNbParts(ind) == maxNb) {
3086         nbFound++;
3087       } else {
3088       }
3089     }
3090     if (nbFound == 0) {
3091       SetErrorCode("There are no blocks, containing the given parts");
3092       return NULL;
3093     }
3094
3095     // Map shapes
3096     TopTools_IndexedMapOfShape anIndices;
3097     TopExp::MapShapes(aBlockOrComp, anIndices);
3098     Handle(TColStd_HArray1OfInteger) anArray;
3099
3100     for (ind = 1; ind <= nbSolids; ind++) {
3101       if (aNbParts(ind) == maxNb) {
3102         anArray = new TColStd_HArray1OfInteger(1,1);
3103         anArray->SetValue(1, anIndices.FindIndex(aSolids(ind)));
3104         anObj = GetEngine()->AddSubShape(theCompound, anArray);
3105         aBlocks->Append(anObj);
3106
3107         // For Python command
3108         TDF_Tool::Entry(anObj->GetEntry(), anEntry);
3109         anAsciiList += anEntry + ", ";
3110         if (aFunction.IsNull())
3111           aFunction = anObj->GetLastFunction();
3112       }
3113     }
3114   }
3115   catch (Standard_Failure) {
3116     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3117     SetErrorCode(aFail->GetMessageString());
3118     return NULL;
3119   }
3120
3121   //Make a Python command
3122   aPartsDescr.Trunc(aPartsDescr.Length() - 2);
3123   anAsciiList.Trunc(anAsciiList.Length() - 2);
3124
3125   GEOM::TPythonDump(aFunction) << "[" << anAsciiList.ToCString()
3126     << "] = geompy.GetBlocksByParts(" << theCompound
3127       << ", [" << aPartsDescr.ToCString() << "])";
3128
3129   SetErrorCode(OK);
3130   return aBlocks;
3131 }
3132
3133 //=============================================================================
3134 /*!
3135  *  MakeMultiTransformation1D
3136  */
3137 //=============================================================================
3138 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation1D
3139                                              (Handle(GEOM_Object)    theObject,
3140                                               const Standard_Integer theDirFace1,
3141                                               const Standard_Integer theDirFace2,
3142                                               const Standard_Integer theNbTimes)
3143 {
3144   SetErrorCode(KO);
3145
3146   if (theObject.IsNull()) return NULL;
3147
3148   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
3149   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be moved
3150
3151   //Add a new Copy object
3152   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
3153
3154   //Add a translate function
3155   Handle(GEOM_Function) aFunction =
3156     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_MULTI_TRANSFORM_1D);
3157
3158   //Check if the function is set correctly
3159   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
3160
3161   GEOMImpl_IBlockTrsf aTI (aFunction);
3162   aTI.SetOriginal(aLastFunction);
3163   aTI.SetFace1U(theDirFace1);
3164   aTI.SetFace2U(theDirFace2);
3165   aTI.SetNbIterU(theNbTimes);
3166
3167   //Compute the transformation
3168   try {
3169 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
3170     OCC_CATCH_SIGNALS;
3171 #endif
3172     if (!GetSolver()->ComputeFunction(aFunction)) {
3173       SetErrorCode("Block driver failed to make multi-transformation");
3174       return NULL;
3175     }
3176   }
3177   catch (Standard_Failure) {
3178     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3179     SetErrorCode(aFail->GetMessageString());
3180     return NULL;
3181   }
3182
3183   //Make a Python command
3184   GEOM::TPythonDump(aFunction) << aCopy << " = geompy.MakeMultiTransformation1D("
3185     << theObject << ", " << theDirFace1 << ", " << theDirFace2 << ", " << theNbTimes << ")";
3186
3187   SetErrorCode(OK);
3188   return aCopy;
3189 }
3190
3191 //=============================================================================
3192 /*!
3193  *  MakeMultiTransformation2D
3194  */
3195 //=============================================================================
3196 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation2D
3197                                              (Handle(GEOM_Object)    theObject,
3198                                               const Standard_Integer theDirFace1U,
3199                                               const Standard_Integer theDirFace2U,
3200                                               const Standard_Integer theNbTimesU,
3201                                               const Standard_Integer theDirFace1V,
3202                                               const Standard_Integer theDirFace2V,
3203                                               const Standard_Integer theNbTimesV)
3204 {
3205   SetErrorCode(KO);
3206
3207   if (theObject.IsNull()) return NULL;
3208
3209   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
3210   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be moved
3211
3212   //Add a new Copy object
3213   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
3214
3215   //Add a translate function
3216   Handle(GEOM_Function) aFunction =
3217     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_MULTI_TRANSFORM_2D);
3218
3219   //Check if the function is set correctly
3220   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
3221
3222   GEOMImpl_IBlockTrsf aTI (aFunction);
3223   aTI.SetOriginal(aLastFunction);
3224   aTI.SetFace1U(theDirFace1U);
3225   aTI.SetFace2U(theDirFace2U);
3226   aTI.SetNbIterU(theNbTimesU);
3227   aTI.SetFace1V(theDirFace1V);
3228   aTI.SetFace2V(theDirFace2V);
3229   aTI.SetNbIterV(theNbTimesV);
3230
3231   //Compute the transformation
3232   try {
3233 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
3234     OCC_CATCH_SIGNALS;
3235 #endif
3236     if (!GetSolver()->ComputeFunction(aFunction)) {
3237       SetErrorCode("Block driver failed to make multi-transformation");
3238       return NULL;
3239     }
3240   }
3241   catch (Standard_Failure) {
3242     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3243     SetErrorCode(aFail->GetMessageString());
3244     return NULL;
3245   }
3246
3247   //Make a Python command
3248   GEOM::TPythonDump(aFunction) << aCopy << " = geompy.MakeMultiTransformation2D("
3249     << theObject << ", " << theDirFace1U << ", " << theDirFace2U << ", " << theNbTimesU
3250       << ", " << theDirFace1V << ", " << theDirFace2V << ", " << theNbTimesV << ")";
3251
3252   SetErrorCode(OK);
3253   return aCopy;
3254 }
3255
3256 //=============================================================================
3257 /*!
3258  *  Propagate
3259  */
3260 //=============================================================================
3261 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::Propagate
3262                                                  (Handle(GEOM_Object) theShape)
3263 {
3264   SetErrorCode(KO);
3265
3266   if (theShape.IsNull()) return NULL;
3267
3268   TopoDS_Shape aShape = theShape->GetValue();
3269   if (aShape.IsNull()) return NULL;
3270
3271   TopTools_IndexedMapOfShape anIndices;
3272   TopExp::MapShapes(aShape, anIndices);
3273
3274   TopTools_IndexedDataMapOfShapeListOfShape MEW;
3275   GEOMImpl_Block6Explorer::MapShapesAndAncestors
3276     (aShape, TopAbs_EDGE, TopAbs_WIRE, MEW);
3277   Standard_Integer ie, nbEdges = MEW.Extent();
3278
3279   // Result
3280   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
3281
3282   TopTools_MapOfShape mapAcceptedEdges;
3283   TCollection_AsciiString aListRes, anEntry;
3284
3285   for (ie = 1; ie <= nbEdges; ie++) {
3286     TopoDS_Shape curE = MEW.FindKey(ie);
3287
3288     if (mapAcceptedEdges.Contains(curE)) continue;
3289
3290     // Build the chain
3291     TopTools_ListOfShape currentChain;
3292     TopTools_ListOfShape listPrevEdges;
3293
3294     currentChain.Append(curE);
3295     listPrevEdges.Append(curE);
3296     mapAcceptedEdges.Add(curE);
3297
3298     // Collect all edges pass by pass
3299     while (listPrevEdges.Extent() > 0) {
3300       // List of edges, added to chain on this cycle pass
3301       TopTools_ListOfShape listCurEdges;
3302
3303       // Find the next portion of edges
3304       TopTools_ListIteratorOfListOfShape itE (listPrevEdges);
3305       for (; itE.More(); itE.Next()) {
3306         TopoDS_Shape anE = itE.Value();
3307
3308         // Iterate on faces, having edge <anE>
3309         TopTools_ListIteratorOfListOfShape itW (MEW.FindFromKey(anE));
3310         for (; itW.More(); itW.Next()) {
3311           TopoDS_Shape aW = itW.Value();
3312           TopoDS_Shape anOppE;
3313
3314           BRepTools_WireExplorer aWE (TopoDS::Wire(aW));
3315           Standard_Integer nb = 1, found = 0;
3316           TopTools_Array1OfShape anEdges (1,4);
3317           for (; aWE.More(); aWE.Next(), nb++) {
3318             if (nb > 4) {
3319               found = 0;
3320               break;
3321             }
3322             anEdges(nb) = aWE.Current();
3323             if (anEdges(nb).IsSame(anE)) found = nb;
3324           }
3325
3326           if (nb == 5 && found > 0) {
3327             // Quadrangle face found, get an opposite edge
3328             Standard_Integer opp = found + 2;
3329             if (opp > 4) opp -= 4;
3330             anOppE = anEdges(opp);
3331
3332             if (!mapAcceptedEdges.Contains(anOppE)) {
3333               // Add found edge to the chain
3334               currentChain.Append(anOppE);
3335               listCurEdges.Append(anOppE);
3336               mapAcceptedEdges.Add(anOppE);
3337             }
3338           } // if (nb == 5 && found > 0)
3339         } // for (; itF.More(); itF.Next())
3340       } // for (; itE.More(); itE.Next())
3341
3342       listPrevEdges = listCurEdges;
3343     } // while (listPrevEdges.Extent() > 0)
3344
3345     // Store the chain in the document
3346     Handle(TColStd_HArray1OfInteger) anArray =
3347       new TColStd_HArray1OfInteger (1, currentChain.Extent());
3348
3349     // Fill array of sub-shape indices
3350     TopTools_ListIteratorOfListOfShape itSub (currentChain);
3351     for (int index = 1; itSub.More(); itSub.Next(), ++index) {
3352       int id = anIndices.FindIndex(itSub.Value());
3353       anArray->SetValue(index, id);
3354     }
3355
3356     // Add a new group object
3357     Handle(GEOM_Object) aChain = GetEngine()->AddSubShape(theShape, anArray);
3358
3359     // Set a GROUP type
3360     aChain->SetType(GEOM_GROUP);
3361
3362     // Set a sub shape type
3363     TDF_Label aFreeLabel = aChain->GetFreeLabel();
3364     TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)TopAbs_EDGE);
3365
3366     // Add the chain to the result
3367     aSeq->Append(aChain);
3368
3369     //Make a Python command
3370     TDF_Tool::Entry(aChain->GetEntry(), anEntry);
3371     aListRes += anEntry + ", ";
3372   }
3373
3374   if (aSeq->IsEmpty()) {
3375     SetErrorCode("There are no quadrangle faces in the shape");
3376     return aSeq;
3377   }
3378
3379   aListRes.Trunc(aListRes.Length() - 2);
3380
3381   // The Propagation doesn't change object so no new function is required.
3382   Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
3383
3384   // Make a Python command
3385   GEOM::TPythonDump(aFunction, /*append=*/true)
3386     << "[" << aListRes.ToCString() << "] = geompy.Propagate(" << theShape << ")";
3387
3388   SetErrorCode(OK);
3389   return aSeq;
3390 }