src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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