src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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