src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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