src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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