src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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