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_Vertex aVert = TopoDS::Vertex(anArg);
 798     TopoDS_Shape aShape = GEOMUtils::GetEdgeNearPoint(aBlockOrComp, aVert);
 799
 800     TopTools_IndexedMapOfShape anIndices;
 801     TopExp::MapShapes(aBlockOrComp, anIndices);
 802     Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
 803     anArray->SetValue(1, anIndices.FindIndex(aShape));
 804     aResult = GetEngine()->AddSubShape(theShape, anArray);
 805   }
 806   catch (Standard_Failure) {
 807     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 808     SetErrorCode(aFail->GetMessageString());
 809     return NULL;
 810   }
 811
 812   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
 813
 814   //Make a Python command
 815   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetEdgeNearPoint("
 816                                << theShape << ", " << thePoint << ")";
 817
 818   SetErrorCode(OK);
 819   return aResult;
 820 }
 821
 822 //=============================================================================
 823 /*!
 824  *  GetFaceByPoints
 825  */
 826 //=============================================================================
 827 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByPoints
 828                                                 (Handle(GEOM_Object) theShape,
 829                                                  Handle(GEOM_Object) thePoint1,
 830                                                  Handle(GEOM_Object) thePoint2,
 831                                                  Handle(GEOM_Object) thePoint3,
 832                                                  Handle(GEOM_Object) thePoint4)
 833 {
 834   SetErrorCode(KO);
 835
 836   //New object
 837   Handle(GEOM_Object) aResult;
 838
 839   // Arguments
 840   if (theShape.IsNull() ||
 841       thePoint1.IsNull() || thePoint2.IsNull() ||
 842       thePoint3.IsNull() || thePoint4.IsNull()) return NULL;
 843
 844   TopoDS_Shape aBlockOrComp = theShape->GetValue();
 845   if (aBlockOrComp.IsNull()) {
 846     SetErrorCode("Block or compound is null");
 847     return NULL;
 848   }
 849
 850   TopoDS_Shape anArg1 = thePoint1->GetValue();
 851   TopoDS_Shape anArg2 = thePoint2->GetValue();
 852   TopoDS_Shape anArg3 = thePoint3->GetValue();
 853   TopoDS_Shape anArg4 = thePoint4->GetValue();
 854   if (anArg1.IsNull() || anArg2.IsNull() ||
 855       anArg3.IsNull() || anArg4.IsNull()) {
 856     SetErrorCode("Null shape is given as argument");
 857     return NULL;
 858   }
 859   if (anArg1.ShapeType() != TopAbs_VERTEX ||
 860       anArg2.ShapeType() != TopAbs_VERTEX ||
 861       anArg3.ShapeType() != TopAbs_VERTEX ||
 862       anArg4.ShapeType() != TopAbs_VERTEX) {
 863     SetErrorCode("Element for face identification is not a vertex");
 864     return NULL;
 865   }
 866
 867   //Compute the Face value
 868   try {
 869 #if OCC_VERSION_LARGE > 0x06010000
 870     OCC_CATCH_SIGNALS;
 871 #endif
 872     TopoDS_Shape aShape;
 873
 874     TopTools_IndexedDataMapOfShapeListOfShape MVF;
 875     GEOMImpl_Block6Explorer::MapShapesAndAncestors(aBlockOrComp, TopAbs_VERTEX, TopAbs_FACE, MVF);
 876
 877     TopoDS_Shape V1,V2,V3,V4;
 878     Standard_Integer ish, ext = MVF.Extent();
 879
 880     if (MVF.Contains(anArg1)) {
 881       V1 = anArg1;
 882     } else {
 883       for (ish = 1; ish <= ext; ish++) {
 884         TopoDS_Shape aShi = MVF.FindKey(ish);
 885         if (BRepTools::Compare(TopoDS::Vertex(anArg1), TopoDS::Vertex(aShi))) {
 886           V1 = aShi;
 887           break;
 888         }
 889       }
 890     }
 891
 892     if (MVF.Contains(anArg2)) {
 893       V2 = anArg2;
 894     } else {
 895       for (ish = 1; ish <= ext; ish++) {
 896         TopoDS_Shape aShi = MVF.FindKey(ish);
 897         if (BRepTools::Compare(TopoDS::Vertex(anArg2), TopoDS::Vertex(aShi))) {
 898           V2 = aShi;
 899           break;
 900         }
 901       }
 902     }
 903
 904     if (MVF.Contains(anArg3)) {
 905       V3 = anArg3;
 906     } else {
 907       for (ish = 1; ish <= ext; ish++) {
 908         TopoDS_Shape aShi = MVF.FindKey(ish);
 909         if (BRepTools::Compare(TopoDS::Vertex(anArg3), TopoDS::Vertex(aShi))) {
 910           V3 = aShi;
 911           break;
 912         }
 913       }
 914     }
 915
 916     if (MVF.Contains(anArg4)) {
 917       V4 = anArg4;
 918     } else {
 919       for (ish = 1; ish <= ext; ish++) {
 920         TopoDS_Shape aShi = MVF.FindKey(ish);
 921         if (BRepTools::Compare(TopoDS::Vertex(anArg4), TopoDS::Vertex(aShi))) {
 922           V4 = aShi;
 923           break;
 924         }
 925       }
 926     }
 927
 928     if (V1.IsNull() || V2.IsNull() || V3.IsNull() || V4.IsNull()) {
 929       SetErrorCode("The given vertex does not belong to the shape");
 930       return NULL;
 931     }
 932
 933     Standard_Integer isFound =
 934       GEOMImpl_Block6Explorer::FindFace(aShape, V1, V2, V3, V4, MVF, Standard_True);
 935     if (isFound == 0) {
 936       SetErrorCode("The given vertices do not belong to one face of the given shape");
 937       return NULL;
 938     } else if (isFound > 1) {
 939       SetErrorCode("The given vertices belong to several faces of the given shape");
 940       return NULL;
 941     } else {
 942       TopTools_IndexedMapOfShape anIndices;
 943       TopExp::MapShapes(aBlockOrComp, anIndices);
 944       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
 945       anArray->SetValue(1, anIndices.FindIndex(aShape));
 946       aResult = GetEngine()->AddSubShape(theShape, anArray);
 947     }
 948   }
 949   catch (Standard_Failure) {
 950     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 951     SetErrorCode(aFail->GetMessageString());
 952     return NULL;
 953   }
 954
 955   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
 956
 957   //Make a Python command
 958   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceByPoints("
 959     << theShape << ", " << thePoint1 << ", " << thePoint2
 960       << ", " << thePoint3 << ", " << thePoint4 << ")";
 961
 962   SetErrorCode(OK);
 963   return aResult;
 964 }
 965
 966 //=============================================================================
 967 /*!
 968  *  GetFaceByEdges
 969  */
 970 //=============================================================================
 971 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByEdges
 972                                                 (Handle(GEOM_Object) theShape,
 973                                                  Handle(GEOM_Object) theEdge1,
 974                                                  Handle(GEOM_Object) theEdge2)
 975 {
 976   SetErrorCode(KO);
 977
 978   //New object
 979   Handle(GEOM_Object) aResult;
 980
 981   // Arguments
 982   if (theShape.IsNull() || theEdge1.IsNull() || theEdge2.IsNull()) return NULL;
 983
 984   TopoDS_Shape aBlockOrComp = theShape->GetValue();
 985   if (aBlockOrComp.IsNull()) {
 986     SetErrorCode("Block or compound is null");
 987     return NULL;
 988   }
 989
 990   TopoDS_Shape anArg1 = theEdge1->GetValue();
 991   TopoDS_Shape anArg2 = theEdge2->GetValue();
 992   if (anArg1.IsNull() || anArg2.IsNull()) {
 993     SetErrorCode("Null shape is given as argument");
 994     return NULL;
 995   }
 996   if (anArg1.ShapeType() != TopAbs_EDGE ||
 997       anArg2.ShapeType() != TopAbs_EDGE) {
 998     SetErrorCode("Element for face identification is not an edge");
 999     return NULL;
1000   }
1001
1002   //Compute the Face value
1003   try {
1004 #if OCC_VERSION_LARGE > 0x06010000
1005     OCC_CATCH_SIGNALS;
1006 #endif
1007     TopoDS_Shape aShape;
1008
1009     TopTools_IndexedDataMapOfShapeListOfShape MEF;
1010     GEOMImpl_Block6Explorer::MapShapesAndAncestors(aBlockOrComp, TopAbs_EDGE, TopAbs_FACE, MEF);
1011
1012     TopoDS_Shape E1,E2;
1013     Standard_Integer ish, ext = MEF.Extent();
1014
1015     if (MEF.Contains(anArg1)) {
1016       E1 = anArg1;
1017     } else {
1018       for (ish = 1; ish <= ext; ish++) {
1019         TopoDS_Shape aShi = MEF.FindKey(ish);
1020         if (GEOMImpl_Block6Explorer::IsSimilarEdges(anArg1, aShi)) {
1021           E1 = aShi;
1022         }
1023       }
1024     }
1025
1026     if (MEF.Contains(anArg2)) {
1027       E2 = anArg2;
1028     } else {
1029       for (ish = 1; ish <= ext; ish++) {
1030         TopoDS_Shape aShi = MEF.FindKey(ish);
1031         if (GEOMImpl_Block6Explorer::IsSimilarEdges(anArg2, aShi)) {
1032           E2 = aShi;
1033         }
1034       }
1035     }
1036
1037     if (E1.IsNull() || E2.IsNull()) {
1038       SetErrorCode("The given edge does not belong to the shape");
1039       return NULL;
1040     }
1041
1042     const TopTools_ListOfShape& aFacesOfE1 = MEF.FindFromKey(E1);
1043     const TopTools_ListOfShape& aFacesOfE2 = MEF.FindFromKey(E2);
1044
1045     Standard_Integer isFound = 0;
1046     TopTools_ListIteratorOfListOfShape anIterF1 (aFacesOfE1);
1047     for (; anIterF1.More(); anIterF1.Next()) {
1048
1049       TopTools_ListIteratorOfListOfShape anIterF2 (aFacesOfE2);
1050       for (; anIterF2.More(); anIterF2.Next()) {
1051
1052         if (anIterF1.Value().IsSame(anIterF2.Value())) {
1053           isFound++;
1054
1055           // Store the face, defined by two edges
1056           aShape = anIterF1.Value();
1057         }
1058       }
1059     }
1060     if (isFound == 0) {
1061       SetErrorCode("The given edges do not belong to one face of the given shape");
1062       return NULL;
1063     } else if (isFound > 1) {
1064       SetErrorCode("The given edges belong to several faces of the given shape");
1065       return NULL;
1066     } else {
1067       TopTools_IndexedMapOfShape anIndices;
1068       TopExp::MapShapes(aBlockOrComp, anIndices);
1069       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
1070       anArray->SetValue(1, anIndices.FindIndex(aShape));
1071       aResult = GetEngine()->AddSubShape(theShape, anArray);
1072     }
1073   }
1074   catch (Standard_Failure) {
1075     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1076     SetErrorCode(aFail->GetMessageString());
1077     return NULL;
1078   }
1079
1080   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1081
1082   //Make a Python command
1083   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceByEdges("
1084     << theShape << ", " << theEdge1 << ", " << theEdge2 << ")";
1085
1086   SetErrorCode(OK);
1087   return aResult;
1088 }
1089
1090 //=============================================================================
1091 /*!
1092  *  GetOppositeFace
1093  */
1094 //=============================================================================
1095 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetOppositeFace
1096                                                 (Handle(GEOM_Object) theShape,
1097                                                  Handle(GEOM_Object) theFace)
1098 {
1099   SetErrorCode(KO);
1100
1101   //New object
1102   Handle(GEOM_Object) aResult;
1103
1104   // Arguments
1105   if (theShape.IsNull() || theFace.IsNull()) return NULL;
1106
1107   TopoDS_Shape aBlockOrComp = theShape->GetValue();
1108   if (aBlockOrComp.IsNull()) {
1109     SetErrorCode("Block is null");
1110     return NULL;
1111   }
1112   if (aBlockOrComp.ShapeType() != TopAbs_SOLID) {
1113     SetErrorCode("Shape is not a block");
1114     return NULL;
1115   }
1116
1117   TopoDS_Shape anArg = theFace->GetValue();
1118   if (anArg.IsNull()) {
1119     SetErrorCode("Null shape is given as argument");
1120     return NULL;
1121   }
1122   if (anArg.ShapeType() != TopAbs_FACE) {
1123     SetErrorCode("Element for face identification is not a face");
1124     return NULL;
1125   }
1126
1127   //Compute the Face value
1128   try {
1129 #if OCC_VERSION_LARGE > 0x06010000
1130     OCC_CATCH_SIGNALS;
1131 #endif
1132     TopoDS_Shape aShape;
1133
1134     GEOMImpl_Block6Explorer aBlockTool;
1135     aBlockTool.InitByBlockAndFace(aBlockOrComp, anArg);
1136     aShape = aBlockTool.GetFace(2);
1137
1138     TopTools_IndexedMapOfShape anIndices;
1139     TopExp::MapShapes(aBlockOrComp, anIndices);
1140     Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
1141     anArray->SetValue(1, anIndices.FindIndex(aShape));
1142     aResult = GetEngine()->AddSubShape(theShape, anArray);
1143   }
1144   catch (Standard_Failure) {
1145     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1146     SetErrorCode(aFail->GetMessageString());
1147     return NULL;
1148   }
1149
1150   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1151
1152   //Make a Python command
1153   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetOppositeFace("
1154                                << theShape << ", " << theFace << ")";
1155
1156   SetErrorCode(OK);
1157   return aResult;
1158 }
1159
1160 //=============================================================================
1161 /*!
1162  *  GetFaceNearPoint
1163  */
1164 //=============================================================================
1165 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceNearPoint
1166                                                 (Handle(GEOM_Object) theShape,
1167                                                  Handle(GEOM_Object) thePoint)
1168 {
1169   SetErrorCode(KO);
1170
1171   //New object
1172   Handle(GEOM_Object) aResult;
1173
1174   // Arguments
1175   if (theShape.IsNull() || thePoint.IsNull()) return NULL;
1176
1177   TopoDS_Shape aBlockOrComp = theShape->GetValue();
1178   if (aBlockOrComp.IsNull()) {
1179     SetErrorCode("Block or compound is null");
1180     return NULL;
1181   }
1182
1183   TopoDS_Shape anArg = thePoint->GetValue();
1184   if (anArg.IsNull()) {
1185     SetErrorCode("Null shape is given as argument");
1186     return NULL;
1187   }
1188   if (anArg.ShapeType() != TopAbs_VERTEX) {
1189     SetErrorCode("Element for face identification is not a vertex");
1190     return NULL;
1191   }
1192
1193   //Compute the Face value
1194   try {
1195 #if OCC_VERSION_LARGE > 0x06010000
1196     OCC_CATCH_SIGNALS;
1197 #endif
1198     TopoDS_Shape aShape;
1199
1200     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
1201     gp_Pnt aPnt = BRep_Tool::Pnt(aVert);
1202     Standard_Real PX, PY, PZ;
1203     aPnt.Coord(PX, PY, PZ);
1204
1205     // 1. Classify the point relatively each face
1206     Standard_Integer nearest = 2, nbFound = 0;
1207     TopTools_DataMapOfShapeInteger mapShapeDist;
1208     TopExp_Explorer exp (aBlockOrComp, TopAbs_FACE);
1209     for (; exp.More(); exp.Next()) {
1210       TopoDS_Shape aFace = exp.Current();
1211
1212       if (!mapShapeDist.IsBound(aFace)) {
1213         Standard_Integer aDistance = 2;
1214
1215         // 1.a. Classify relatively Surface
1216         Handle(Geom_Surface) aSurf = BRep_Tool::Surface(TopoDS::Face(aFace));
1217         Handle(ShapeAnalysis_Surface) aSurfAna = new ShapeAnalysis_Surface (aSurf);
1218         gp_Pnt2d p2dOnSurf = aSurfAna->ValueOfUV(aPnt, Precision::Confusion());
1219         gp_Pnt p3dOnSurf = aSurfAna->Value(p2dOnSurf);
1220         Standard_Real aDist = p3dOnSurf.Distance(aPnt);
1221         if (aDist > Precision::Confusion()) {
1222           // OUT of Surface
1223           aDistance = 1;
1224         } else {
1225           // 1.b. Classify relatively the face itself
1226           BRepClass_FaceClassifier FC (TopoDS::Face(aFace), p2dOnSurf, Precision::Confusion());
1227           if (FC.State() == TopAbs_IN) {
1228             aDistance = -1;
1229           } else if (FC.State() == TopAbs_ON) {
1230             aDistance = 0;
1231           } else { // OUT
1232             aDistance = 1;
1233           }
1234         }
1235
1236         if (aDistance < nearest) {
1237           nearest = aDistance;
1238           aShape = aFace;
1239           nbFound = 1;
1240
1241           // A first found face, containing the point inside, will be returned.
1242           // It is the solution, if there are no
1243           // coincident or intersecting faces in the compound.
1244           if (nearest == -1) break;
1245
1246         } else if (aDistance == nearest) {
1247           nbFound++;
1248         } else {
1249         }
1250
1251         mapShapeDist.Bind(aFace, aDistance);
1252       } // if (!mapShapeDist.IsBound(aFace))
1253     }
1254
1255     // 2. Define face, containing the point or having minimum distance to it
1256     if (nbFound > 1) {
1257       if (nearest == 0) {
1258         // The point is on boundary of some faces and there are
1259         // no faces, having the point inside
1260         SetErrorCode("Multiple faces near the given point are found");
1261         return NULL;
1262
1263       } else if (nearest == 1) {
1264         // The point is outside some faces and there are
1265         // no faces, having the point inside or on boundary.
1266         // We will get a nearest face
1267         Standard_Real bigReal = RealLast();
1268         Standard_Real minDist = bigReal;
1269         TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
1270         for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
1271           if (mapShapeDistIter.Value() == 1) {
1272             TopoDS_Shape aFace = mapShapeDistIter.Key();
1273             Standard_Real aDist = bigReal;
1274
1275             // 2.a. Fast check of distance - if point projection on surface is on face
1276             Handle(Geom_Surface) aSurf = BRep_Tool::Surface(TopoDS::Face(aFace));
1277             Handle(ShapeAnalysis_Surface) aSurfAna = new ShapeAnalysis_Surface (aSurf);
1278             gp_Pnt2d p2dOnSurf = aSurfAna->ValueOfUV(aPnt, Precision::Confusion());
1279             gp_Pnt p3dOnSurf = aSurfAna->Value(p2dOnSurf);
1280             aDist = p3dOnSurf.Distance(aPnt);
1281
1282             BRepClass_FaceClassifier FC (TopoDS::Face(aFace), p2dOnSurf, Precision::Confusion());
1283             if (FC.State() == TopAbs_OUT) {
1284               if (aDist < minDist) {
1285                 // 2.b. Slow check - if point projection on surface is outside of face
1286                 BRepExtrema_DistShapeShape aDistTool (aVert, aFace);
1287                 if (!aDistTool.IsDone()) {
1288                   SetErrorCode("Can not find a distance from the given point to one of faces");
1289                   return NULL;
1290                 }
1291                 aDist = aDistTool.Value();
1292               } else {
1293                 aDist = bigReal;
1294               }
1295             }
1296
1297             if (aDist < minDist) {
1298               minDist = aDist;
1299               aShape = aFace;
1300             }
1301           }
1302         }
1303       } else { // nearest == -1
1304 //        // The point is inside some faces.
1305 //        // We will get a face with nearest center
1306 //        Standard_Real minDist = RealLast();
1307 //        TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
1308 //        for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
1309 //          if (mapShapeDistIter.Value() == -1) {
1310 //            TopoDS_Shape aFace = mapShapeDistIter.Key();
1311 //            GProp_GProps aSystem;
1312 //            BRepGProp::SurfaceProperties(aFace, aSystem);
1313 //            gp_Pnt aCenterMass = aSystem.CentreOfMass();
1314 //
1315 //            Standard_Real aDist = aCenterMass.Distance(aPnt);
1316 //            if (aDist < minDist) {
1317 //              minDist = aDist;
1318 //              aShape = aFace;
1319 //            }
1320 //          }
1321 //        }
1322       }
1323     } // if (nbFound > 1)
1324
1325     if (nbFound == 0) {
1326       SetErrorCode("There are no faces near the given point");
1327       return NULL;
1328     } else {
1329       TopTools_IndexedMapOfShape anIndices;
1330       TopExp::MapShapes(aBlockOrComp, anIndices);
1331       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
1332       anArray->SetValue(1, anIndices.FindIndex(aShape));
1333       aResult = GetEngine()->AddSubShape(theShape, anArray);
1334     }
1335   }
1336   catch (Standard_Failure) {
1337     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1338     SetErrorCode(aFail->GetMessageString());
1339     return NULL;
1340   }
1341
1342   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1343
1344   //Make a Python command
1345   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceNearPoint("
1346                                << theShape << ", " << thePoint << ")";
1347
1348   SetErrorCode(OK);
1349   return aResult;
1350 }
1351
1352 //=============================================================================
1353 /*!
1354  *  GetFaceByNormale
1355  */
1356 //=============================================================================
1357 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByNormale
1358                                                 (Handle(GEOM_Object) theShape,
1359                                                  Handle(GEOM_Object) theVector)
1360 {
1361   SetErrorCode(KO);
1362
1363   //New object
1364   Handle(GEOM_Object) aResult;
1365
1366   // Arguments
1367   if (theShape.IsNull() || theVector.IsNull()) return NULL;
1368
1369   TopoDS_Shape aBlockOrComp = theShape->GetValue();
1370   if (aBlockOrComp.IsNull()) {
1371     SetErrorCode("Block or compound is null");
1372     return NULL;
1373   }
1374
1375   TopoDS_Shape anArg = theVector->GetValue();
1376   if (anArg.IsNull()) {
1377     SetErrorCode("Null shape is given as argument");
1378     return NULL;
1379   }
1380   if (anArg.ShapeType() != TopAbs_EDGE) {
1381     SetErrorCode("Element for normale identification is not an edge");
1382     return NULL;
1383   }
1384
1385   //Compute the Face value
1386   try {
1387 #if OCC_VERSION_LARGE > 0x06010000
1388     OCC_CATCH_SIGNALS;
1389 #endif
1390     TopoDS_Shape aShape;
1391
1392     TopoDS_Edge anEdge = TopoDS::Edge(anArg);
1393     TopoDS_Vertex V1, V2;
1394     TopExp::Vertices(anEdge, V1, V2, Standard_True);
1395     gp_Pnt P1 = BRep_Tool::Pnt(V1);
1396     gp_Pnt P2 = BRep_Tool::Pnt(V2);
1397     gp_Vec aVec (P1, P2);
1398     if (aVec.Magnitude() < Precision::Confusion()) {
1399       SetErrorCode("Vector with null magnitude is given");
1400       return NULL;
1401     }
1402
1403     Standard_Real minAngle = RealLast();
1404     TopTools_MapOfShape mapShape;
1405     TopExp_Explorer exp (aBlockOrComp, TopAbs_FACE);
1406     for (; exp.More(); exp.Next()) {
1407       if (mapShape.Add(exp.Current())) {
1408         TopoDS_Face aFace = TopoDS::Face(exp.Current());
1409         BRepAdaptor_Surface SF (aFace);
1410
1411         Standard_Real u, v, x;
1412
1413         // find a point on the surface to get normal direction in
1414         u = SF.FirstUParameter();
1415         x = SF.LastUParameter();
1416         if (Precision::IsInfinite(u)) {
1417           u =  (Precision::IsInfinite(x)) ? 0. : x;
1418         } else if (!Precision::IsInfinite(x)) {
1419           u = (u+x) / 2.;
1420         }
1421
1422         v = SF.FirstVParameter();
1423         x = SF.LastVParameter();
1424         if (Precision::IsInfinite(v)) {
1425           v =  (Precision::IsInfinite(x)) ? 0. : x;
1426         } else if (!Precision::IsInfinite(x)) {
1427           v = (v+x) / 2.;
1428         }
1429
1430         // compute the normal direction
1431         gp_Vec Vec1,Vec2;
1432         SF.D1(u,v,P1,Vec1,Vec2);
1433         gp_Vec V = Vec1.Crossed(Vec2);
1434         x = V.Magnitude();
1435         if (V.Magnitude() < Precision::Confusion()) {
1436           SetErrorCode("Normal vector of a face has null magnitude");
1437           return NULL;
1438         }
1439
1440         // consider the face orientation
1441         if (aFace.Orientation() == TopAbs_REVERSED ||
1442             aFace.Orientation() == TopAbs_INTERNAL) {
1443           V = - V;
1444         }
1445
1446         // compute the angle and compare with the minimal one
1447         Standard_Real anAngle = aVec.Angle(V);
1448         if (anAngle < minAngle) {
1449           minAngle = anAngle;
1450           aShape = aFace;
1451         }
1452       }
1453     }
1454
1455     if (aShape.IsNull()) {
1456       SetErrorCode("Failed to find a face by the given normale");
1457       return NULL;
1458     } else {
1459       TopTools_IndexedMapOfShape anIndices;
1460       TopExp::MapShapes(aBlockOrComp, anIndices);
1461       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
1462       anArray->SetValue(1, anIndices.FindIndex(aShape));
1463       aResult = GetEngine()->AddSubShape(theShape, anArray);
1464     }
1465   }
1466   catch (Standard_Failure) {
1467     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1468     SetErrorCode(aFail->GetMessageString());
1469     return NULL;
1470   }
1471
1472   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1473
1474   //Make a Python command
1475   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceByNormale("
1476     << theShape << ", " << theVector << ")";
1477
1478   SetErrorCode(OK);
1479   return aResult;
1480 }
1481
1482 //=============================================================================
1483 /*!
1484  *  GetShapesNearPoint
1485  */
1486 //=============================================================================
1487 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetShapesNearPoint
1488                                          (Handle(GEOM_Object)    theShape,
1489                                           Handle(GEOM_Object)    thePoint,
1490                                           const Standard_Integer theShapeType,
1491                                           const Standard_Real    theConstTolerance)
1492 {
1493   SetErrorCode(KO);
1494
1495   // New object
1496   Handle(GEOM_Object) aResult;
1497
1498   // Arguments
1499   if (theShape.IsNull() || thePoint.IsNull()) return NULL;
1500
1501   TopoDS_Shape aBlockOrComp = theShape->GetValue();
1502   if (aBlockOrComp.IsNull()) {
1503     SetErrorCode("Block or compound is null");
1504     return NULL;
1505   }
1506
1507   TopoDS_Shape anArg = thePoint->GetValue();
1508   if (anArg.IsNull()) {
1509     SetErrorCode("Null shape is given as argument");
1510     return NULL;
1511   }
1512   if (anArg.ShapeType() != TopAbs_VERTEX) {
1513     SetErrorCode("Element for face identification is not a vertex");
1514     return NULL;
1515   }
1516
1517   if (theShapeType < TopAbs_SOLID || TopAbs_VERTEX < theShapeType) {
1518     SetErrorCode("Invalid type of result is requested");
1519     return NULL;
1520   }
1521
1522   Standard_Real theTolerance = theConstTolerance;
1523   if (theTolerance < Precision::Confusion()) {
1524     theTolerance = Precision::Confusion();
1525   }
1526
1527   // Compute the result
1528   try {
1529 #if OCC_VERSION_LARGE > 0x06010000
1530     OCC_CATCH_SIGNALS;
1531 #endif
1532     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
1533
1534     TopTools_MapOfShape mapShape;
1535     Standard_Integer nbEdges = 0;
1536     TopExp_Explorer exp (aBlockOrComp, TopAbs_ShapeEnum(theShapeType));
1537     for (; exp.More(); exp.Next()) {
1538       if (mapShape.Add(exp.Current())) {
1539         nbEdges++;
1540       }
1541     }
1542
1543     if (nbEdges == 0) {
1544       SetErrorCode("Given shape contains no sub-shapes of requested type");
1545       return NULL;
1546     }
1547
1548     // Calculate distances and find min
1549     mapShape.Clear();
1550     Standard_Integer ind = 1;
1551     Standard_Real aMinDist = RealLast();
1552     TopTools_Array1OfShape anEdges (1, nbEdges);
1553     TColStd_Array1OfReal aDistances (1, nbEdges);
1554     for (exp.Init(aBlockOrComp, TopAbs_ShapeEnum(theShapeType)); exp.More(); exp.Next()) {
1555       if (mapShape.Add(exp.Current())) {
1556         TopoDS_Shape anEdge = exp.Current();
1557         anEdges(ind) = anEdge;
1558
1559         BRepExtrema_DistShapeShape aDistTool (aVert, anEdges(ind));
1560         if (!aDistTool.IsDone()) {
1561           SetErrorCode("Can not find a distance from the given point to one of sub-shapes");
1562           return NULL;
1563         }
1564         aDistances(ind) = aDistTool.Value();
1565         if (aDistances(ind) < aMinDist) {
1566           aMinDist = aDistances(ind);
1567         }
1568         ind++;
1569       }
1570     }
1571
1572     if (aMinDist < RealLast()) {
1573       // Collect sub-shapes with distance < (aMinDist + theTolerance)
1574       int nbSubShapes = 0;
1575       TopTools_Array1OfShape aNearShapes (1, nbEdges);
1576       for (ind = 1; ind <= nbEdges; ind++) {
1577         if (aDistances(ind) < aMinDist + theTolerance) {
1578           nbSubShapes++;
1579           aNearShapes(nbSubShapes) = anEdges(ind);
1580         }
1581       }
1582
1583       // Add sub-shape
1584       TopTools_IndexedMapOfShape anIndices;
1585       TopExp::MapShapes(aBlockOrComp, anIndices);
1586       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger (1, nbSubShapes);
1587       for (ind = 1; ind <= nbSubShapes; ind++) {
1588         anArray->SetValue(ind, anIndices.FindIndex(aNearShapes(ind)));
1589       }
1590       aResult = GetEngine()->AddSubShape(theShape, anArray);
1591     }
1592   }
1593   catch (Standard_Failure) {
1594     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1595     SetErrorCode(aFail->GetMessageString());
1596     return NULL;
1597   }
1598
1599   if (aResult.IsNull())
1600     return NULL;
1601
1602   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1603
1604   //Make a Python command
1605   GEOM::TPythonDump(aFunction)
1606     << aResult << " = geompy.GetShapesNearPoint(" << theShape << ", " << thePoint
1607     << ", " << TopAbs_ShapeEnum(theShapeType) << ", " << theTolerance << ")";
1608
1609   SetErrorCode(OK);
1610   return aResult;
1611 }
1612
1613 //=============================================================================
1614 /*!
1615  *  IsCompoundOfBlocks
1616  */
1617 //=============================================================================
1618 Standard_Boolean GEOMImpl_IBlocksOperations::IsCompoundOfBlocks
1619                                                 (Handle(GEOM_Object)    theCompound,
1620                                                  const Standard_Integer theMinNbFaces,
1621                                                  const Standard_Integer theMaxNbFaces,
1622                                                  Standard_Integer&      theNbBlocks)
1623 {
1624   SetErrorCode(KO);
1625   Standard_Boolean isCompOfBlocks = Standard_False;
1626   theNbBlocks = 0;
1627
1628   if (theCompound.IsNull()) return isCompOfBlocks;
1629   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
1630
1631   //Check
1632   isCompOfBlocks = Standard_True;
1633   try {
1634 #if OCC_VERSION_LARGE > 0x06010000
1635     OCC_CATCH_SIGNALS;
1636 #endif
1637     TopTools_MapOfShape mapShape;
1638     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
1639     for (; exp.More(); exp.Next()) {
1640       if (mapShape.Add(exp.Current())) {
1641         TopoDS_Shape aSolid = exp.Current();
1642
1643         TopTools_MapOfShape mapFaces;
1644         TopExp_Explorer expF (aSolid, TopAbs_FACE);
1645         Standard_Integer nbFaces = 0;
1646         for (; expF.More(); expF.Next()) {
1647           if (mapFaces.Add(expF.Current())) {
1648             nbFaces++;
1649             if (nbFaces > theMaxNbFaces) {
1650               isCompOfBlocks = Standard_False;
1651               break;
1652             }
1653           }
1654         }
1655         if (nbFaces < theMinNbFaces || theMaxNbFaces < nbFaces) {
1656           isCompOfBlocks = Standard_False;
1657         } else {
1658           theNbBlocks++;
1659         }
1660       }
1661     }
1662   }
1663   catch (Standard_Failure) {
1664     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1665     SetErrorCode(aFail->GetMessageString());
1666     return isCompOfBlocks;
1667   }
1668
1669   SetErrorCode(OK);
1670   return isCompOfBlocks;
1671 }
1672
1673 //=============================================================================
1674 /*!
1675  *  Set of functions, used by CheckCompoundOfBlocks() method
1676  */
1677 //=============================================================================
1678 void GEOMImpl_IBlocksOperations::AddBlocksFrom (const TopoDS_Shape&   theShape,
1679                                                 TopTools_ListOfShape& BLO,
1680                                                 TopTools_ListOfShape& NOT,
1681                                                 TopTools_ListOfShape& EXT,
1682                                                 TopTools_ListOfShape& NOQ)
1683 {
1684   TopAbs_ShapeEnum aType = theShape.ShapeType();
1685   switch (aType) {
1686   case TopAbs_COMPOUND:
1687   case TopAbs_COMPSOLID:
1688     {
1689       TopoDS_Iterator It (theShape);
1690       for (; It.More(); It.Next()) {
1691         AddBlocksFrom(It.Value(), BLO, NOT, EXT, NOQ);
1692       }
1693     }
1694     break;
1695   case TopAbs_SOLID:
1696     {
1697       // Check, if there are seam or degenerated edges
1698       BlockFix_CheckTool aTool;
1699       aTool.SetShape(theShape);
1700       aTool.Perform();
1701       if (aTool.NbPossibleBlocks() > 0) {
1702         EXT.Append(theShape);
1703       } else {
1704         // Count faces and edges in each face to recognize blocks
1705         TopTools_MapOfShape mapFaces;
1706         Standard_Integer nbFaces = 0;
1707         Standard_Boolean hasNonQuadr = Standard_False;
1708         TopExp_Explorer expF (theShape, TopAbs_FACE);
1709
1710         for (; expF.More(); expF.Next()) {
1711           if (mapFaces.Add(expF.Current())) {
1712             nbFaces++;
1713             //0021483//if (nbFaces > 6) break;
1714
1715             // get wire
1716             TopoDS_Shape aF = expF.Current();
1717             TopExp_Explorer wires (aF, TopAbs_WIRE);
1718             if (!wires.More()) {
1719               // no wire in the face
1720               hasNonQuadr = Standard_True;
1721               NOQ.Append(aF);//0021483
1722               //0021483//break;
1723               continue;
1724             }
1725             TopoDS_Shape aWire = wires.Current();
1726             wires.Next();
1727             if (wires.More()) {
1728               // multiple wires in the face
1729               hasNonQuadr = Standard_True;
1730               NOQ.Append(aF);//0021483
1731               //0021483//break;
1732               continue;
1733             }
1734
1735             // Check number of edges in the face
1736             Standard_Integer nbEdges = 0;
1737             TopTools_MapOfShape mapEdges;
1738             TopExp_Explorer expW (aWire, TopAbs_EDGE);
1739             for (; expW.More(); expW.Next()) {
1740               if (mapEdges.Add(expW.Current())) {
1741                 nbEdges++;
1742                 if (nbEdges > 4) break;
1743               }
1744             }
1745             if (nbEdges != 4) {
1746               hasNonQuadr = Standard_True;
1747               NOQ.Append(aF);//0021483
1748             }
1749           }
1750         }
1751
1752         if (nbFaces == 6 && !hasNonQuadr) {
1753           BLO.Append(theShape);
1754         } else {
1755           NOT.Append(theShape);
1756         }
1757       }
1758     }
1759     break;
1760   case TopAbs_SHELL: //0021483
1761   case TopAbs_FACE: //0021483
1762     {
1763       // Count edges in each face
1764       TopTools_MapOfShape mapFaces;
1765       TopExp_Explorer expF (theShape, TopAbs_FACE);
1766       for (; expF.More(); expF.Next()) {
1767         if (mapFaces.Add(expF.Current())) {
1768           // get wire
1769           TopoDS_Shape aF = expF.Current();
1770           TopExp_Explorer wires (aF, TopAbs_WIRE);
1771           if (!wires.More()) {
1772             // no wire in the face
1773             NOQ.Append(aF);//0021483
1774             continue;
1775           }
1776           TopoDS_Shape aWire = wires.Current();
1777           wires.Next();
1778           if (wires.More()) {
1779             // multiple wires in the face
1780             NOQ.Append(aF);//0021483
1781             continue;
1782           }
1783
1784           // Check number of edges in the face
1785           Standard_Integer nbEdges = 0;
1786           TopTools_MapOfShape mapEdges;
1787           TopExp_Explorer expW (aWire, TopAbs_EDGE);
1788           for (; expW.More(); expW.Next()) {
1789             if (mapEdges.Add(expW.Current())) {
1790               nbEdges++;
1791               if (nbEdges > 4) break;
1792             }
1793           }
1794           if (nbEdges != 4) {
1795             NOQ.Append(aF);//0021483
1796           }
1797         }
1798       }
1799     }
1800     break;
1801   default:
1802     NOT.Append(theShape);
1803   }
1804 }
1805
1806 void AddBlocksFromOld (const TopoDS_Shape&   theShape,
1807                        TopTools_ListOfShape& BLO,
1808                        TopTools_ListOfShape& NOT,
1809                        TopTools_ListOfShape& DEG,
1810                        TopTools_ListOfShape& SEA)
1811 {
1812   TopAbs_ShapeEnum aType = theShape.ShapeType();
1813   switch (aType) {
1814   case TopAbs_COMPOUND:
1815   case TopAbs_COMPSOLID:
1816     {
1817       TopoDS_Iterator It (theShape);
1818       for (; It.More(); It.Next()) {
1819         AddBlocksFromOld(It.Value(), BLO, NOT, DEG, SEA);
1820       }
1821     }
1822     break;
1823   case TopAbs_SOLID:
1824     {
1825       TopTools_MapOfShape mapFaces;
1826       TopExp_Explorer expF (theShape, TopAbs_FACE);
1827       Standard_Integer nbFaces = 0;
1828       Standard_Boolean hasNonQuadr = Standard_False;
1829       Standard_Boolean hasDegenerated = Standard_False;
1830       Standard_Boolean hasSeam = Standard_False;
1831       for (; expF.More(); expF.Next()) {
1832         if (mapFaces.Add(expF.Current())) {
1833           nbFaces++;
1834           if (nbFaces > 6) break;
1835
1836           // Check number of edges in the face
1837           Standard_Integer nbEdges = 0;
1838           TopTools_MapOfShape mapEdges;
1839
1840           // get wire
1841           TopoDS_Shape aF = expF.Current();
1842           TopExp_Explorer wires (aF, TopAbs_WIRE);
1843           if (!wires.More()) {
1844             // no wire in the face
1845             hasNonQuadr = Standard_True;
1846             break;
1847           }
1848           TopoDS_Shape aWire = wires.Current();
1849           wires.Next();
1850           if (wires.More()) {
1851             // multiple wires in the face
1852             hasNonQuadr = Standard_True;
1853             break;
1854           }
1855
1856           // iterate on wire
1857           BRepTools_WireExplorer aWE (TopoDS::Wire(aWire), TopoDS::Face(aF));
1858           for (; aWE.More(); aWE.Next(), nbEdges++) {
1859             if (BRep_Tool::Degenerated(aWE.Current())) {
1860               // degenerated edge found
1861               hasDegenerated = Standard_True;
1862 //              break;
1863             }
1864             if (mapEdges.Contains(aWE.Current())) {
1865               // seam edge found
1866               hasSeam = Standard_True;
1867 //              break;
1868             }
1869             mapEdges.Add(aWE.Current());
1870           }
1871           if (nbEdges != 4) {
1872             hasNonQuadr = Standard_True;
1873           }
1874         }
1875       }
1876       if (nbFaces == 6) {
1877         if (hasDegenerated || hasSeam) {
1878           if (hasDegenerated) {
1879             DEG.Append(theShape);
1880           }
1881           if (hasSeam) {
1882             SEA.Append(theShape);
1883           }
1884         } else if (hasNonQuadr) {
1885           NOT.Append(theShape);
1886         } else {
1887           BLO.Append(theShape);
1888         }
1889       } else {
1890         NOT.Append(theShape);
1891       }
1892     }
1893     break;
1894   default:
1895     NOT.Append(theShape);
1896   }
1897 }
1898
1899 #define REL_NOT_CONNECTED 0
1900 #define REL_OK            1
1901 #define REL_NOT_GLUED     2
1902 #define REL_COLLISION_VV  3
1903 #define REL_COLLISION_FF  4
1904 #define REL_COLLISION_EE  5
1905 #define REL_UNKNOWN       6
1906
1907 Standard_Integer BlocksRelation (const TopoDS_Shape& theBlock1,
1908                                  const TopoDS_Shape& theBlock2)
1909 {
1910   // Compare bounding boxes before calling BRepExtrema_DistShapeShape
1911   Standard_Real Xmin1, Ymin1, Zmin1, Xmax1, Ymax1, Zmax1;
1912   Standard_Real Xmin2, Ymin2, Zmin2, Xmax2, Ymax2, Zmax2;
1913   Bnd_Box B1, B2;
1914   BRepBndLib::Add(theBlock1, B1);
1915   BRepBndLib::Add(theBlock2, B2);
1916   B1.Get(Xmin1, Ymin1, Zmin1, Xmax1, Ymax1, Zmax1);
1917   B2.Get(Xmin2, Ymin2, Zmin2, Xmax2, Ymax2, Zmax2);
1918   if (Xmax2 < Xmin1 || Xmax1 < Xmin2 ||
1919       Ymax2 < Ymin1 || Ymax1 < Ymin2 ||
1920       Zmax2 < Zmin1 || Zmax1 < Zmin2) {
1921     return REL_NOT_CONNECTED;
1922   }
1923
1924   BRepExtrema_DistShapeShape dst (theBlock1, theBlock2);
1925   if (!dst.IsDone()) {
1926     return REL_UNKNOWN;
1927   }
1928
1929   if (dst.Value() > Precision::Confusion()) {
1930     return REL_NOT_CONNECTED;
1931   }
1932
1933   if (dst.InnerSolution()) {
1934     return REL_COLLISION_VV;
1935   }
1936
1937   Standard_Integer nbSol = dst.NbSolution();
1938   Standard_Integer relation = REL_OK;
1939   Standard_Integer nbVerts = 0;
1940   Standard_Integer nbEdges = 0;
1941   Standard_Integer sol = 1;
1942   for (; sol <= nbSol; sol++) {
1943     BRepExtrema_SupportType supp1 = dst.SupportTypeShape1(sol);
1944     BRepExtrema_SupportType supp2 = dst.SupportTypeShape2(sol);
1945     if (supp1 == BRepExtrema_IsVertex && supp2 == BRepExtrema_IsVertex) {
1946       nbVerts++;
1947     } else if (supp1 == BRepExtrema_IsInFace || supp2 == BRepExtrema_IsInFace) {
1948       return REL_COLLISION_FF;
1949     } else if (supp1 == BRepExtrema_IsOnEdge && supp2 == BRepExtrema_IsOnEdge) {
1950       nbEdges++;
1951     } else if ((supp1 == BRepExtrema_IsOnEdge && supp2 == BRepExtrema_IsVertex) ||
1952                (supp2 == BRepExtrema_IsOnEdge && supp1 == BRepExtrema_IsVertex)) {
1953       relation = REL_COLLISION_EE;
1954     } else {
1955     }
1956   }
1957
1958   if (relation != REL_OK) {
1959     return relation;
1960   }
1961
1962   TColStd_Array1OfInteger vertSol (1, nbVerts);
1963   TopTools_Array1OfShape V1 (1, nbVerts);
1964   TopTools_Array1OfShape V2 (1, nbVerts);
1965   Standard_Integer ivs = 0;
1966   for (sol = 1; sol <= nbSol; sol++) {
1967     if (dst.SupportTypeShape1(sol) == BRepExtrema_IsVertex &&
1968         dst.SupportTypeShape2(sol) == BRepExtrema_IsVertex) {
1969       TopoDS_Vertex Vcur = TopoDS::Vertex(dst.SupportOnShape1(sol));
1970       // Check, that this vertex is far enough from other solution vertices.
1971       Standard_Integer ii = 1;
1972       for (; ii <= ivs; ii++) {
1973         if (BRepTools::Compare(TopoDS::Vertex(V1(ii)), Vcur)) {
1974           continue;
1975         }
1976       }
1977       ivs++;
1978       vertSol(ivs) = sol;
1979       V1(ivs) = Vcur;
1980       V2(ivs) = dst.SupportOnShape2(sol);
1981     }
1982   }
1983
1984   // As we deal only with quadrangles,
1985   // 2, 3 or 4 vertex solutions can be found.
1986   if (ivs <= 1) {
1987     if (nbEdges > 0) {
1988       return REL_COLLISION_FF;
1989     }
1990     return REL_NOT_CONNECTED;
1991   }
1992   if (ivs > 4) {
1993     return REL_UNKNOWN;
1994   }
1995
1996   // Check sharing of coincident entities.
1997   if (ivs == 2 || ivs == 3) {
1998     // Map vertices and edges of the blocks
1999     TopTools_IndexedDataMapOfShapeListOfShape MVE1, MVE2;
2000     GEOMImpl_Block6Explorer::MapShapesAndAncestors
2001       (theBlock1, TopAbs_VERTEX, TopAbs_EDGE, MVE1);
2002     GEOMImpl_Block6Explorer::MapShapesAndAncestors
2003       (theBlock2, TopAbs_VERTEX, TopAbs_EDGE, MVE2);
2004
2005     if (ivs == 2) {
2006       // Find common edge
2007       TopoDS_Shape anEdge1, anEdge2;
2008       GEOMImpl_Block6Explorer::FindEdge(anEdge1, V1(1), V1(2), MVE1);
2009       if (anEdge1.IsNull()) return REL_UNKNOWN;
2010
2011       GEOMImpl_Block6Explorer::FindEdge(anEdge2, V2(1), V2(2), MVE2);
2012       if (anEdge2.IsNull()) return REL_UNKNOWN;
2013
2014       if (!anEdge1.IsSame(anEdge2)) return REL_NOT_GLUED;
2015
2016     } else { // ivs == 3
2017       // Find common edges
2018       Standard_Integer e1_v1 = 1;
2019       Standard_Integer e1_v2 = 2;
2020       Standard_Integer e2_v1 = 3;
2021       Standard_Integer e2_v2 = 1;
2022
2023       TopoDS_Shape anEdge11, anEdge12;
2024       GEOMImpl_Block6Explorer::FindEdge(anEdge11, V1(e1_v1), V1(e1_v2), MVE1);
2025       if (anEdge11.IsNull()) {
2026         e1_v2 = 3;
2027         e2_v1 = 2;
2028         GEOMImpl_Block6Explorer::FindEdge(anEdge11, V1(e1_v1), V1(e1_v2), MVE1);
2029         if (anEdge11.IsNull()) return REL_UNKNOWN;
2030       }
2031       GEOMImpl_Block6Explorer::FindEdge(anEdge12, V1(e2_v1), V1(e2_v2), MVE1);
2032       if (anEdge12.IsNull()) {
2033         e2_v2 = 5 - e2_v1;
2034         GEOMImpl_Block6Explorer::FindEdge(anEdge12, V1(e2_v1), V1(e2_v2), MVE1);
2035         if (anEdge12.IsNull()) return REL_UNKNOWN;
2036       }
2037
2038       TopoDS_Shape anEdge21, anEdge22;
2039       GEOMImpl_Block6Explorer::FindEdge(anEdge21, V2(e1_v1), V2(e1_v2), MVE2);
2040       if (anEdge21.IsNull()) return REL_UNKNOWN;
2041       GEOMImpl_Block6Explorer::FindEdge(anEdge22, V2(e2_v1), V2(e2_v2), MVE2);
2042       if (anEdge22.IsNull()) return REL_UNKNOWN;
2043
2044       // Check of edges coincidence (with some precision) have to be done here
2045       // if (!anEdge11.IsEqual(anEdge21)) return REL_UNKNOWN;
2046       // if (!anEdge12.IsEqual(anEdge22)) return REL_UNKNOWN;
2047
2048       // Check of edges sharing
2049       if (!anEdge11.IsSame(anEdge21)) return REL_NOT_GLUED;
2050       if (!anEdge12.IsSame(anEdge22)) return REL_NOT_GLUED;
2051     }
2052   }
2053
2054   if (ivs == 4) {
2055     // Map vertices and faces of the blocks
2056     TopTools_IndexedDataMapOfShapeListOfShape MVF1, MVF2;
2057     GEOMImpl_Block6Explorer::MapShapesAndAncestors
2058       (theBlock1, TopAbs_VERTEX, TopAbs_FACE, MVF1);
2059     GEOMImpl_Block6Explorer::MapShapesAndAncestors
2060       (theBlock2, TopAbs_VERTEX, TopAbs_FACE, MVF2);
2061
2062     TopoDS_Shape aFace1, aFace2;
2063     GEOMImpl_Block6Explorer::FindFace(aFace1, V1(1), V1(2), V1(3), V1(4), MVF1);
2064     if (aFace1.IsNull()) return REL_UNKNOWN;
2065     GEOMImpl_Block6Explorer::FindFace(aFace2, V2(1), V2(2), V2(3), V2(4), MVF2);
2066     if (aFace2.IsNull()) return REL_UNKNOWN;
2067
2068     // Check of faces coincidence (with some precision) have to be done here
2069     // if (!aFace1.IsEqual(aFace2)) return REL_UNKNOWN;
2070
2071     // Check of faces sharing
2072     if (!aFace1.IsSame(aFace2)) return REL_NOT_GLUED;
2073   }
2074
2075   return REL_OK;
2076 }
2077
2078 void FindConnected (const Standard_Integer         theBlockIndex,
2079                     const TColStd_Array2OfInteger& theRelations,
2080                     TColStd_MapOfInteger&          theProcessedMap,
2081                     TColStd_MapOfInteger&          theConnectedMap)
2082 {
2083   theConnectedMap.Add(theBlockIndex);
2084   theProcessedMap.Add(theBlockIndex);
2085
2086   Standard_Integer nbBlocks = theRelations.ColLength();
2087   Standard_Integer col = 1;
2088   for (; col <= nbBlocks; col++) {
2089     if (theRelations(theBlockIndex, col) == REL_OK ||
2090         theRelations(theBlockIndex, col) == REL_NOT_GLUED) {
2091       if (!theProcessedMap.Contains(col)) {
2092         FindConnected(col, theRelations, theProcessedMap, theConnectedMap);
2093       }
2094     }
2095   }
2096 }
2097
2098 Standard_Boolean HasAnyConnection (const Standard_Integer         theBlockIndex,
2099                                    const TColStd_MapOfInteger&    theWith,
2100                                    const TColStd_Array2OfInteger& theRelations,
2101                                    TColStd_MapOfInteger&          theProcessedMap)
2102 {
2103   theProcessedMap.Add(theBlockIndex);
2104
2105   Standard_Integer nbBlocks = theRelations.ColLength();
2106   Standard_Integer col = 1;
2107   for (; col <= nbBlocks; col++) {
2108     if (theRelations(theBlockIndex, col) != REL_NOT_CONNECTED) {
2109       if (!theProcessedMap.Contains(col)) {
2110         if (theWith.Contains(col))
2111           return Standard_True;
2112         if (HasAnyConnection(col, theWith, theRelations, theProcessedMap))
2113           return Standard_True;
2114       }
2115     }
2116   }
2117
2118   return Standard_False;
2119 }
2120
2121 //=============================================================================
2122 /*!
2123  *  CheckCompoundOfBlocksOld
2124  */
2125 //=============================================================================
2126 Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocksOld
2127                                                 (Handle(GEOM_Object) theCompound,
2128                                                  std::list<BCError>&      theErrors)
2129 {
2130   SetErrorCode(KO);
2131
2132   if (theCompound.IsNull()) return Standard_False;
2133   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2134
2135   Standard_Boolean isCompOfBlocks = Standard_True;
2136
2137   // Map sub-shapes and their indices
2138   TopTools_IndexedMapOfShape anIndices;
2139   TopExp::MapShapes(aBlockOrComp, anIndices);
2140
2141   // 1. Report non-blocks
2142   TopTools_ListOfShape NOT; // Not blocks
2143   TopTools_ListOfShape DEG; // Hexahedral solids, having degenerated edges
2144   TopTools_ListOfShape SEA; // Hexahedral solids, having seam edges
2145   TopTools_ListOfShape BLO; // All blocks from the given compound
2146   AddBlocksFromOld(aBlockOrComp, BLO, NOT, DEG, SEA);
2147
2148   if (NOT.Extent() > 0) {
2149     isCompOfBlocks = Standard_False;
2150     BCError anErr;
2151     anErr.error = NOT_BLOCK;
2152     TopTools_ListIteratorOfListOfShape it (NOT);
2153     for (; it.More(); it.Next()) {
2154       anErr.incriminated.push_back(anIndices.FindIndex(it.Value()));
2155     }
2156     theErrors.push_back(anErr);
2157   }
2158
2159   if (DEG.Extent() > 0 || SEA.Extent() > 0) {
2160     isCompOfBlocks = Standard_False;
2161     BCError anErr;
2162     anErr.error = EXTRA_EDGE;
2163
2164     TopTools_ListIteratorOfListOfShape itDEG (DEG);
2165     for (; itDEG.More(); itDEG.Next()) {
2166       anErr.incriminated.push_back(anIndices.FindIndex(itDEG.Value()));
2167     }
2168
2169     TopTools_ListIteratorOfListOfShape itSEA (SEA);
2170     for (; itSEA.More(); itSEA.Next()) {
2171       anErr.incriminated.push_back(anIndices.FindIndex(itSEA.Value()));
2172     }
2173
2174     theErrors.push_back(anErr);
2175   }
2176
2177   Standard_Integer nbBlocks = BLO.Extent();
2178   if (nbBlocks == 0) {
2179     isCompOfBlocks = Standard_False;
2180     SetErrorCode(OK);
2181     return isCompOfBlocks;
2182   }
2183   if (nbBlocks == 1) {
2184     SetErrorCode(OK);
2185     return isCompOfBlocks;
2186   }
2187
2188   // Convert list of blocks into array for easy and fast access
2189   Standard_Integer ibl = 1;
2190   TopTools_Array1OfShape aBlocks (1, nbBlocks);
2191   TopTools_ListIteratorOfListOfShape BLOit (BLO);
2192   for (; BLOit.More(); BLOit.Next(), ibl++) {
2193     aBlocks.SetValue(ibl, BLOit.Value());
2194   }
2195
2196   // 2. Find relations between all blocks,
2197   //    report connection errors (NOT_GLUED and INVALID_CONNECTION)
2198   TColStd_Array2OfInteger aRelations (1, nbBlocks, 1, nbBlocks);
2199   aRelations.Init(REL_NOT_CONNECTED);
2200
2201   Standard_Integer row = 1;
2202   for (row = 1; row <= nbBlocks; row++) {
2203     TopoDS_Shape aBlock = aBlocks.Value(row);
2204
2205     Standard_Integer col = row + 1;
2206     for (; col <= nbBlocks; col++) {
2207       Standard_Integer aRel = BlocksRelation(aBlock, aBlocks.Value(col));
2208       if (aRel != REL_NOT_CONNECTED) {
2209         aRelations.SetValue(row, col, aRel);
2210         aRelations.SetValue(col, row, aRel);
2211         if (aRel == REL_NOT_GLUED) {
2212           // report connection error
2213           isCompOfBlocks = Standard_False;
2214           BCError anErr;
2215           anErr.error = NOT_GLUED;
2216           anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(row)));
2217           anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(col)));
2218           theErrors.push_back(anErr);
2219         } else if (aRel == REL_COLLISION_VV ||
2220                    aRel == REL_COLLISION_FF ||
2221                    aRel == REL_COLLISION_EE ||
2222                    aRel == REL_UNKNOWN) {
2223           // report connection error
2224           isCompOfBlocks = Standard_False;
2225           BCError anErr;
2226           anErr.error = INVALID_CONNECTION;
2227           anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(row)));
2228           anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(col)));
2229           theErrors.push_back(anErr);
2230         } else {
2231         }
2232       }
2233     }
2234   }
2235
2236   // 3. Find largest set of connected (good connection or not glued) blocks
2237   TColStd_MapOfInteger aProcessedMap;
2238   TColStd_MapOfInteger aLargestSet;
2239   TColStd_MapOfInteger aCurrentSet;
2240   for (ibl = 1; ibl <= nbBlocks; ibl++) {
2241     if (!aProcessedMap.Contains(ibl)) {
2242       aCurrentSet.Clear();
2243       FindConnected(ibl, aRelations, aProcessedMap, aCurrentSet);
2244       if (aCurrentSet.Extent() > aLargestSet.Extent()) {
2245         aLargestSet = aCurrentSet;
2246       }
2247     }
2248   }
2249
2250   // 4. Report all blocks, isolated from <aLargestSet>
2251   BCError anErr;
2252   anErr.error = NOT_CONNECTED;
2253   Standard_Boolean hasIsolated = Standard_False;
2254   for (ibl = 1; ibl <= nbBlocks; ibl++) {
2255     if (!aLargestSet.Contains(ibl)) {
2256       aProcessedMap.Clear();
2257       if (!HasAnyConnection(ibl, aLargestSet, aRelations, aProcessedMap)) {
2258         // report connection absence
2259         hasIsolated = Standard_True;
2260         anErr.incriminated.push_back(anIndices.FindIndex(aBlocks.Value(ibl)));
2261       }
2262     }
2263   }
2264   if (hasIsolated) {
2265     isCompOfBlocks = Standard_False;
2266     theErrors.push_back(anErr);
2267   }
2268
2269   SetErrorCode(OK);
2270   return isCompOfBlocks;
2271 }
2272
2273 //=============================================================================
2274 /*!
2275  *  PrintBCErrors
2276  */
2277 //=============================================================================
2278 TCollection_AsciiString GEOMImpl_IBlocksOperations::PrintBCErrors
2279                                                 (Handle(GEOM_Object)  theCompound,
2280                                                  const std::list<BCError>& theErrors)
2281 {
2282   TCollection_AsciiString aDescr;
2283
2284   std::list<BCError>::const_iterator errIt = theErrors.begin();
2285   int i = 0;
2286   for (; errIt != theErrors.end(); i++, errIt++) {
2287     BCError errStruct = *errIt;
2288
2289     switch (errStruct.error) {
2290     case NOT_BLOCK:
2291       aDescr += "\n\tNot a Blocks: ";
2292       break;
2293     case EXTRA_EDGE:
2294       aDescr += "\n\tHexahedral solids with degenerated and/or seam edges: ";
2295       break;
2296     case INVALID_CONNECTION:
2297       aDescr += "\n\tInvalid connection between two blocks: ";
2298       break;
2299     case NOT_CONNECTED:
2300       aDescr += "\n\tBlocks, not connected with main body: ";
2301       break;
2302     case NOT_GLUED:
2303       aDescr += "\n\tNot glued blocks: ";
2304       break;
2305     default:
2306       break;
2307     }
2308
2309     std::list<int> sshList = errStruct.incriminated;
2310     std::list<int>::iterator sshIt = sshList.begin();
2311     int jj = 0;
2312     for (; sshIt != sshList.end(); jj++, sshIt++) {
2313       if (jj > 0)
2314         aDescr += ", ";
2315       aDescr += TCollection_AsciiString(*sshIt);
2316     }
2317   }
2318
2319   return aDescr;
2320 }
2321
2322 //=============================================================================
2323 /*!
2324  *  CheckCompoundOfBlocks
2325  */
2326 //=============================================================================
2327 Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
2328                                               (Handle(GEOM_Object) theCompound,
2329                                                std::list<BCError>& theErrors)
2330 {
2331   SetErrorCode(KO);
2332
2333   if (theCompound.IsNull()) return Standard_False;
2334   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2335
2336   Standard_Boolean isCompOfBlocks = Standard_True;
2337
2338   // Map sub-shapes and their indices
2339   TopTools_IndexedMapOfShape anIndices;
2340   TopExp::MapShapes(aBlockOrComp, anIndices);
2341
2342   // 1. Separate blocks from non-blocks
2343   TopTools_ListOfShape NOT; // Not blocks
2344   TopTools_ListOfShape EXT; // Hexahedral solids, having degenerated and/or seam edges
2345   TopTools_ListOfShape BLO; // All blocks from the given compound
2346   TopTools_ListOfShape NOQ; // All non-quadrangular faces
2347   AddBlocksFrom(aBlockOrComp, BLO, NOT, EXT, NOQ);
2348
2349   // Report non-blocks
2350   if (NOT.Extent() > 0) {
2351     isCompOfBlocks = Standard_False;
2352     BCError anErr;
2353     anErr.error = NOT_BLOCK;
2354     TopTools_ListIteratorOfListOfShape it (NOT);
2355     for (; it.More(); it.Next()) {
2356       anErr.incriminated.push_back(anIndices.FindIndex(it.Value()));
2357     }
2358     theErrors.push_back(anErr);
2359   }
2360
2361   // Report solids, having degenerated and/or seam edges
2362   if (EXT.Extent() > 0) {
2363     isCompOfBlocks = Standard_False;
2364     BCError anErr;
2365     anErr.error = EXTRA_EDGE;
2366     TopTools_ListIteratorOfListOfShape it (EXT);
2367     for (; it.More(); it.Next()) {
2368       anErr.incriminated.push_back(anIndices.FindIndex(it.Value()));
2369     }
2370     theErrors.push_back(anErr);
2371   }
2372
2373   Standard_Integer nbBlocks = BLO.Extent();
2374   if (nbBlocks == 0) {
2375     isCompOfBlocks = Standard_False;
2376     SetErrorCode(OK);
2377     return isCompOfBlocks;
2378   }
2379   if (nbBlocks == 1) {
2380     SetErrorCode(OK);
2381     return isCompOfBlocks;
2382   }
2383
2384   // Prepare data for 2. and 3.
2385   TColStd_Array2OfInteger aRelations (1, nbBlocks, 1, nbBlocks);
2386   aRelations.Init(REL_NOT_CONNECTED);
2387
2388   TopTools_IndexedMapOfShape mapBlocks;
2389
2390   BRep_Builder BB;
2391   TopoDS_Compound aComp;
2392   BB.MakeCompound(aComp);
2393
2394   TopTools_ListIteratorOfListOfShape BLOit (BLO);
2395   for (; BLOit.More(); BLOit.Next()) {
2396     mapBlocks.Add(BLOit.Value());
2397     BB.Add(aComp, BLOit.Value());
2398   }
2399
2400   // 2. Find glued blocks (having shared faces)
2401   TopTools_IndexedDataMapOfShapeListOfShape mapFaceBlocks;
2402   GEOMImpl_Block6Explorer::MapShapesAndAncestors
2403     (aComp, TopAbs_FACE, TopAbs_SOLID, mapFaceBlocks);
2404
2405   Standard_Integer prevInd = 0, curInd = 0;
2406   Standard_Integer ind = 1, nbFaces = mapFaceBlocks.Extent();
2407   for (; ind <= nbFaces; ind++) {
2408     const TopTools_ListOfShape& aGluedBlocks = mapFaceBlocks.FindFromIndex(ind);
2409     if (aGluedBlocks.Extent() > 1) { // Shared face found
2410       TopTools_ListIteratorOfListOfShape aGluedBlocksIt (aGluedBlocks);
2411       TopoDS_Shape prevBlock, curBlock;
2412       for (; aGluedBlocksIt.More(); aGluedBlocksIt.Next()) {
2413         curBlock = aGluedBlocksIt.Value();
2414         if (!prevBlock.IsNull()) {
2415           prevInd = mapBlocks.FindIndex(prevBlock);
2416           curInd  = mapBlocks.FindIndex(curBlock);
2417           aRelations.SetValue(prevInd, curInd, REL_OK);
2418           aRelations.SetValue(curInd, prevInd, REL_OK);
2419         }
2420         prevBlock = curBlock;
2421       }
2422     }
2423   }
2424
2425   // 3. Find not glued blocks
2426   GEOMAlgo_GlueAnalyser aGD;
2427
2428   aGD.SetShape(aComp);
2429   aGD.SetTolerance(Precision::Confusion());
2430   aGD.SetCheckGeometry(Standard_True);
2431   aGD.Perform();
2432
2433   Standard_Integer iErr, iWrn;
2434   iErr = aGD.ErrorStatus();
2435   if (iErr) {
2436     SetErrorCode("Error in GEOMAlgo_GlueAnalyser");
2437     return isCompOfBlocks;
2438   }
2439   iWrn = aGD.WarningStatus();
2440   if (iWrn) {
2441     MESSAGE("Warning in GEOMAlgo_GlueAnalyser");
2442   }
2443
2444   // Report not glued blocks
2445   if (aGD.HasSolidsToGlue()) {
2446     isCompOfBlocks = Standard_False;
2447     Standard_Integer aSx1Ind, aSx2Ind;
2448
2449     const GEOMAlgo_ListOfCoupleOfShapes& aLCS = aGD.SolidsToGlue();
2450     GEOMAlgo_ListIteratorOfListOfCoupleOfShapes aItCS (aLCS);
2451     for (; aItCS.More(); aItCS.Next()) {
2452       const GEOMAlgo_CoupleOfShapes& aCS = aItCS.Value();
2453       const TopoDS_Shape& aSx1 = aCS.Shape1();
2454       const TopoDS_Shape& aSx2 = aCS.Shape2();
2455
2456       aSx1Ind = mapBlocks.FindIndex(aSx1);
2457       aSx2Ind = mapBlocks.FindIndex(aSx2);
2458       aRelations.SetValue(aSx1Ind, aSx2Ind, NOT_GLUED);
2459       aRelations.SetValue(aSx2Ind, aSx1Ind, NOT_GLUED);
2460
2461       BCError anErr;
2462       anErr.error = NOT_GLUED;
2463       anErr.incriminated.push_back(anIndices.FindIndex(aSx1));
2464       anErr.incriminated.push_back(anIndices.FindIndex(aSx2));
2465       theErrors.push_back(anErr);
2466     }
2467   }
2468
2469   // 4. Find largest set of connected (good connection or not glued) blocks
2470   Standard_Integer ibl = 1;
2471   TColStd_MapOfInteger aProcessedMap;
2472   TColStd_MapOfInteger aLargestSet;
2473   TColStd_MapOfInteger aCurrentSet;
2474   for (ibl = 1; ibl <= nbBlocks; ibl++) {
2475     if (!aProcessedMap.Contains(ibl)) {
2476       aCurrentSet.Clear();
2477       FindConnected(ibl, aRelations, aProcessedMap, aCurrentSet);
2478       if (aCurrentSet.Extent() > aLargestSet.Extent()) {
2479         aLargestSet = aCurrentSet;
2480       }
2481     }
2482   }
2483
2484   // 5. Report all blocks, isolated from <aLargestSet>
2485   BCError anErr;
2486   anErr.error = NOT_CONNECTED;
2487   Standard_Boolean hasIsolated = Standard_False;
2488   for (ibl = 1; ibl <= nbBlocks; ibl++) {
2489     if (!aLargestSet.Contains(ibl)) {
2490       aProcessedMap.Clear();
2491       if (!HasAnyConnection(ibl, aLargestSet, aRelations, aProcessedMap)) {
2492         // report connection absence
2493         hasIsolated = Standard_True;
2494         anErr.incriminated.push_back(anIndices.FindIndex(mapBlocks.FindKey(ibl)));
2495       }
2496     }
2497   }
2498   if (hasIsolated) {
2499     isCompOfBlocks = Standard_False;
2500     theErrors.push_back(anErr);
2501   }
2502
2503   SetErrorCode(OK);
2504   return isCompOfBlocks;
2505 }
2506
2507 //=============================================================================
2508 /*!
2509  *  GetNonBlocks
2510  */
2511 //=============================================================================
2512 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetNonBlocks
2513                                      (Handle(GEOM_Object) theShape,
2514                                       Handle(GEOM_Object)& theNonQuads)
2515 {
2516   SetErrorCode(KO);
2517
2518   if (theShape.IsNull()) return NULL;
2519   TopoDS_Shape aShape = theShape->GetValue();
2520
2521   // Separate blocks from non-blocks
2522   TopTools_ListOfShape BLO; // All blocks from the given compound
2523   TopTools_ListOfShape NOT; // Not blocks
2524   TopTools_ListOfShape EXT; // Hexahedral solids, having degenerated and/or seam edges
2525   TopTools_ListOfShape NOQ; // All non-quadrangular faces
2526   AddBlocksFrom(aShape, BLO, NOT, EXT, NOQ);
2527
2528   if (NOT.IsEmpty() && EXT.IsEmpty() && NOQ.IsEmpty()) {
2529     SetErrorCode("NOT_FOUND_ANY");
2530     return NULL;
2531   }
2532
2533   // Map sub-shapes and their indices
2534   TopTools_IndexedMapOfShape anIndices;
2535   TopExp::MapShapes(aShape, anIndices);
2536
2537   // Non-blocks
2538   Handle(GEOM_Object) aNonBlocks;
2539   if (NOT.Extent() > 0 || EXT.Extent() > 0) {
2540     Handle(TColStd_HArray1OfInteger) anArray =
2541       new TColStd_HArray1OfInteger (1, NOT.Extent() + EXT.Extent());
2542     Standard_Integer ii = 1;
2543     TopTools_ListIteratorOfListOfShape it1 (NOT);
2544     for (; it1.More(); it1.Next(), ii++) {
2545       anArray->SetValue(ii, anIndices.FindIndex(it1.Value()));
2546     }
2547     TopTools_ListIteratorOfListOfShape it2 (EXT);
2548     for (; it2.More(); it2.Next(), ii++) {
2549       anArray->SetValue(ii, anIndices.FindIndex(it2.Value()));
2550     }
2551
2552     aNonBlocks = GetEngine()->AddSubShape(theShape, anArray);
2553     if (aNonBlocks.IsNull()) {
2554       SetErrorCode("Error in algorithm: result found, but cannot be returned.");
2555       return NULL;
2556     }
2557     aNonBlocks->SetType(GEOM_GROUP);
2558     TDF_Label aFreeLabel = aNonBlocks->GetFreeLabel();
2559     TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)TopAbs_SOLID);
2560   }
2561
2562   // Non-quadrangles
2563   if (NOQ.Extent() > 0) {
2564     Handle(TColStd_HArray1OfInteger) anArray =
2565       new TColStd_HArray1OfInteger (1, NOQ.Extent());
2566     Standard_Integer ii = 1;
2567     TopTools_ListIteratorOfListOfShape it1 (NOQ);
2568     for (; it1.More(); it1.Next(), ii++) {
2569       anArray->SetValue(ii, anIndices.FindIndex(it1.Value()));
2570     }
2571
2572     theNonQuads = GetEngine()->AddSubShape(theShape, anArray);
2573     if (theNonQuads.IsNull()) {
2574       SetErrorCode("Error in algorithm: result found, but cannot be returned.");
2575       return NULL;
2576     }
2577     theNonQuads->SetType(GEOM_GROUP);
2578     TDF_Label aFreeLabel = theNonQuads->GetFreeLabel();
2579     TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)TopAbs_FACE);
2580   }
2581
2582   //Make a Python command
2583   Handle(GEOM_Function) aMainShape = theShape->GetLastFunction();
2584   GEOM::TPythonDump pd (aMainShape, /*append=*/true);
2585   pd << "(";
2586   if (aNonBlocks.IsNull())
2587     pd << "no_bad_solids";
2588   else
2589     pd << aNonBlocks;
2590   pd << ", ";
2591   if (theNonQuads.IsNull())
2592     pd << "no_bad_faces";
2593   else
2594     pd << theNonQuads;
2595   pd << ") = geompy.GetNonBlocks(" << theShape << ")";
2596
2597   SetErrorCode(OK);
2598   return aNonBlocks;
2599 }
2600
2601 //=============================================================================
2602 /*!
2603  *  RemoveExtraEdges
2604  */
2605 //=============================================================================
2606 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::RemoveExtraEdges
2607                                      (Handle(GEOM_Object) theObject,
2608                                       const Standard_Integer theOptimumNbFaces)
2609 {
2610   SetErrorCode(KO);
2611
2612   if (theObject.IsNull()) return NULL;
2613
2614   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
2615   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be fixed
2616
2617   //Add a new Copy object
2618   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
2619
2620   //Add a function
2621   Handle(GEOM_Function) aFunction =
2622     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_REMOVE_EXTRA);
2623
2624   //Check if the function is set correctly
2625   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
2626
2627   GEOMImpl_IBlockTrsf aTI (aFunction);
2628   aTI.SetOriginal(aLastFunction);
2629   aTI.SetOptimumNbFaces(theOptimumNbFaces);
2630
2631   //Compute the fixed shape
2632   try {
2633 #if OCC_VERSION_LARGE > 0x06010000
2634     OCC_CATCH_SIGNALS;
2635 #endif
2636     if (!GetSolver()->ComputeFunction(aFunction)) {
2637       SetErrorCode("Block driver failed to remove extra edges of the given shape");
2638       return NULL;
2639     }
2640   }
2641   catch (Standard_Failure) {
2642     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2643     SetErrorCode(aFail->GetMessageString());
2644     return NULL;
2645   }
2646
2647   //Make a Python command
2648   std::string doUnionFaces = (theOptimumNbFaces < 0) ? "False" : "True";
2649   GEOM::TPythonDump(aFunction) << aCopy << " = geompy.RemoveExtraEdges("
2650                                << theObject << ", " << doUnionFaces.data() << ")";
2651
2652   SetErrorCode(OK);
2653   return aCopy;
2654 }
2655
2656 //=============================================================================
2657 /*!
2658  *  CheckAndImprove
2659  */
2660 //=============================================================================
2661 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::CheckAndImprove
2662                                              (Handle(GEOM_Object) theObject)
2663 {
2664   SetErrorCode(KO);
2665
2666   if (theObject.IsNull()) return NULL;
2667
2668   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
2669   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be fixed
2670
2671   //Add a new Copy object
2672   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
2673
2674   //Add a function
2675   Handle(GEOM_Function) aFunction =
2676     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_COMPOUND_IMPROVE);
2677
2678   //Check if the function is set correctly
2679   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
2680
2681   GEOMImpl_IBlockTrsf aTI (aFunction);
2682   aTI.SetOriginal(aLastFunction);
2683
2684   // -1 means do not unite faces on common surface (?except case of seam edge between them?)
2685   //aTI.SetOptimumNbFaces(-1);
2686   aTI.SetOptimumNbFaces(6);
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 improve the given blocks compound");
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   GEOM::TPythonDump(aFunction) << aCopy
2706     << " = geompy.CheckAndImprove(" << theObject << ")";
2707
2708   SetErrorCode(OK);
2709   return aCopy;
2710 }
2711
2712 //=============================================================================
2713 /*!
2714  *  ExplodeCompoundOfBlocks
2715  */
2716 //=============================================================================
2717 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::ExplodeCompoundOfBlocks
2718                                                 (Handle(GEOM_Object)    theCompound,
2719                                                  const Standard_Integer theMinNbFaces,
2720                                                  const Standard_Integer theMaxNbFaces)
2721 {
2722   SetErrorCode(KO);
2723
2724   if (theCompound.IsNull()) return NULL;
2725   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2726   if (aBlockOrComp.IsNull()) return NULL;
2727
2728   Handle(TColStd_HSequenceOfTransient) aBlocks = new TColStd_HSequenceOfTransient;
2729   Handle(GEOM_Object) anObj;
2730   Handle(GEOM_Function) aFunction;
2731
2732   TopTools_MapOfShape mapShape;
2733   TCollection_AsciiString anAsciiList, anEntry;
2734
2735   // Map shapes
2736   TopTools_IndexedMapOfShape anIndices;
2737   TopExp::MapShapes(aBlockOrComp, anIndices);
2738   Handle(TColStd_HArray1OfInteger) anArray;
2739
2740   // Explode
2741   try {
2742 #if OCC_VERSION_LARGE > 0x06010000
2743     OCC_CATCH_SIGNALS;
2744 #endif
2745     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
2746     for (; exp.More(); exp.Next()) {
2747       if (mapShape.Add(exp.Current())) {
2748         TopoDS_Shape aSolid = exp.Current();
2749
2750         TopTools_MapOfShape mapFaces;
2751         TopExp_Explorer expF (aSolid, TopAbs_FACE);
2752         Standard_Integer nbFaces = 0;
2753         for (; expF.More(); expF.Next()) {
2754           if (mapFaces.Add(expF.Current())) {
2755             nbFaces++;
2756           }
2757         }
2758
2759         if (theMinNbFaces <= nbFaces && nbFaces <= theMaxNbFaces) {
2760           anArray = new TColStd_HArray1OfInteger(1,1);
2761           anArray->SetValue(1, anIndices.FindIndex(aSolid));
2762           anObj = GetEngine()->AddSubShape(theCompound, anArray);
2763           aBlocks->Append(anObj);
2764
2765           //Make a Python command
2766           TDF_Tool::Entry(anObj->GetEntry(), anEntry);
2767           anAsciiList += anEntry + ", ";
2768         }
2769       }
2770     }
2771   }
2772   catch (Standard_Failure) {
2773     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2774     SetErrorCode(aFail->GetMessageString());
2775     return aBlocks;
2776   }
2777
2778   if (aBlocks->IsEmpty()) {
2779     SetErrorCode("There are no specified blocks in the given shape");
2780     return aBlocks;
2781   }
2782
2783   anAsciiList.Trunc(anAsciiList.Length() - 2);
2784
2785   //The explode doesn't change object so no new function is required.
2786   aFunction = theCompound->GetLastFunction();
2787
2788   //Make a Python command
2789   GEOM::TPythonDump(aFunction, /*append=*/true)
2790     << "[" << anAsciiList.ToCString() << "] = geompy.MakeBlockExplode("
2791     << theCompound << ", " << theMinNbFaces << ", " << theMaxNbFaces << ")";
2792
2793   SetErrorCode(OK);
2794   return aBlocks;
2795 }
2796
2797 //=============================================================================
2798 /*!
2799  *  GetBlockNearPoint
2800  */
2801 //=============================================================================
2802 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockNearPoint
2803                                                 (Handle(GEOM_Object) theCompound,
2804                                                  Handle(GEOM_Object) thePoint)
2805 {
2806   SetErrorCode(KO);
2807
2808   //New object
2809   Handle(GEOM_Object) aResult;
2810
2811   // Arguments
2812   if (theCompound.IsNull() || thePoint.IsNull()) return NULL;
2813
2814   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2815   if (aBlockOrComp.IsNull()) {
2816     SetErrorCode("Compound is null");
2817     return NULL;
2818   }
2819   if (aBlockOrComp.ShapeType() != TopAbs_COMPOUND &&
2820       aBlockOrComp.ShapeType() != TopAbs_COMPSOLID) {
2821     SetErrorCode("Shape to find block in is not a compound");
2822     return NULL;
2823   }
2824
2825   TopoDS_Shape anArg = thePoint->GetValue();
2826   if (anArg.IsNull()) {
2827     SetErrorCode("Point is null");
2828     return NULL;
2829   }
2830   if (anArg.ShapeType() != TopAbs_VERTEX) {
2831     SetErrorCode("Shape for block identification is not a vertex");
2832     return NULL;
2833   }
2834
2835   //Compute the Block value
2836   try {
2837 #if OCC_VERSION_LARGE > 0x06010000
2838     OCC_CATCH_SIGNALS;
2839 #endif
2840     TopoDS_Shape aShape;
2841
2842     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
2843     gp_Pnt aPnt = BRep_Tool::Pnt(aVert);
2844     Standard_Real PX, PY, PZ;
2845     aPnt.Coord(PX, PY, PZ);
2846
2847     // 1. Classify the point relatively each block
2848     Standard_Integer nearest = 2, nbFound = 0;
2849     TopTools_DataMapOfShapeInteger mapShapeDist;
2850     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
2851     for (; exp.More(); exp.Next()) {
2852       TopoDS_Shape aSolid = exp.Current();
2853
2854       if (!mapShapeDist.IsBound(aSolid)) {
2855         Standard_Integer aDistance = 2;
2856
2857         // 1.a. Classify relatively Bounding box
2858         Standard_Real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
2859         Bnd_Box BB;
2860         BRepBndLib::Add(aSolid, BB);
2861         BB.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
2862         if (PX < Xmin || Xmax < PX ||
2863             PY < Ymin || Ymax < PY ||
2864             PZ < Zmin || Zmax < PZ) {
2865           // OUT of bounding box
2866           aDistance = 1;
2867         } else {
2868           // 1.b. Classify relatively the solid itself
2869           BRepClass3d_SolidClassifier SC (aSolid, aPnt, Precision::Confusion());
2870           if (SC.State() == TopAbs_IN) {
2871             aDistance = -1;
2872           } else if (SC.State() == TopAbs_ON) {
2873             aDistance = 0;
2874           } else { // OUT
2875             aDistance = 1;
2876           }
2877         }
2878
2879         if (aDistance < nearest) {
2880           nearest = aDistance;
2881           aShape = aSolid;
2882           nbFound = 1;
2883
2884           // A first found block, containing the point inside, will be returned.
2885           // It is the solution, if there are no intersecting blocks in the compound.
2886           if (nearest == -1) break;
2887
2888         } else if (aDistance == nearest) {
2889           nbFound++;
2890         } else {
2891         }
2892
2893         mapShapeDist.Bind(aSolid, aDistance);
2894       } // if (!mapShapeDist.IsBound(aSolid))
2895     }
2896
2897     // 2. Define block, containing the point or having minimum distance to it
2898     if (nbFound > 1) {
2899       if (nearest == 0) {
2900         // The point is on boundary of some blocks and there are
2901         // no blocks, having the point inside their volume
2902         SetErrorCode("Multiple blocks near the given point are found");
2903         return NULL;
2904
2905       } else if (nearest == 1) {
2906         // The point is outside some blocks and there are
2907         // no blocks, having the point inside or on boundary.
2908         // We will get a nearest block
2909         Standard_Real minDist = RealLast();
2910         TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
2911         for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
2912           if (mapShapeDistIter.Value() == 1) {
2913             TopoDS_Shape aSolid = mapShapeDistIter.Key();
2914             BRepExtrema_DistShapeShape aDistTool (aVert, aSolid);
2915             if (!aDistTool.IsDone()) {
2916               SetErrorCode("Can not find a distance from the given point to one of blocks");
2917               return NULL;
2918             }
2919             Standard_Real aDist = aDistTool.Value();
2920             if (aDist < minDist) {
2921               minDist = aDist;
2922               aShape = aSolid;
2923             }
2924           }
2925         }
2926       } else { // nearest == -1
2927 //        // The point is inside some blocks.
2928 //        // We will get a block with nearest center
2929 //        Standard_Real minDist = RealLast();
2930 //        TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
2931 //        for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
2932 //          if (mapShapeDistIter.Value() == -1) {
2933 //            TopoDS_Shape aSolid = mapShapeDistIter.Key();
2934 //            GProp_GProps aSystem;
2935 //            BRepGProp::VolumeProperties(aSolid, aSystem);
2936 //            gp_Pnt aCenterMass = aSystem.CentreOfMass();
2937 //
2938 //            Standard_Real aDist = aCenterMass.Distance(aPnt);
2939 //            if (aDist < minDist) {
2940 //              minDist = aDist;
2941 //              aShape = aSolid;
2942 //            }
2943 //          }
2944 //        }
2945       }
2946     } // if (nbFound > 1)
2947
2948     if (nbFound == 0) {
2949       SetErrorCode("There are no blocks near the given point");
2950       return NULL;
2951     } else {
2952       TopTools_IndexedMapOfShape anIndices;
2953       TopExp::MapShapes(aBlockOrComp, anIndices);
2954       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
2955       anArray->SetValue(1, anIndices.FindIndex(aShape));
2956       aResult = GetEngine()->AddSubShape(theCompound, anArray);
2957     }
2958   }
2959   catch (Standard_Failure) {
2960     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2961     SetErrorCode(aFail->GetMessageString());
2962     return NULL;
2963   }
2964
2965   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
2966
2967   //Make a Python command
2968   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetBlockNearPoint("
2969                                << theCompound << ", " << thePoint << ")";
2970
2971   SetErrorCode(OK);
2972   return aResult;
2973 }
2974
2975 //=============================================================================
2976 /*!
2977  *  GetBlockByParts
2978  */
2979 //=============================================================================
2980 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockByParts
2981                       (Handle(GEOM_Object)                         theCompound,
2982                        const Handle(TColStd_HSequenceOfTransient)& theParts)
2983 {
2984   SetErrorCode(KO);
2985
2986   Handle(GEOM_Object) aResult;
2987
2988   if (theCompound.IsNull() || theParts.IsNull()) return NULL;
2989   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2990   if (aBlockOrComp.IsNull()) return NULL;
2991
2992   //Get the parts
2993   Standard_Integer argi, aLen = theParts->Length();
2994   TopTools_Array1OfShape anArgs (1, aLen);
2995   TCollection_AsciiString anEntry, aPartsDescr;
2996   for (argi = 1; argi <= aLen; argi++) {
2997     Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast(theParts->Value(argi));
2998     Handle(GEOM_Function) aRef = anObj->GetLastFunction();
2999     if (aRef.IsNull()) return NULL;
3000
3001     TopoDS_Shape anArg = aRef->GetValue();
3002     if (anArg.IsNull()) {
3003       SetErrorCode("Null shape is given as argument");
3004       return NULL;
3005     }
3006     anArgs(argi) = anArg;
3007
3008     // For Python command
3009     TDF_Tool::Entry(anObj->GetEntry(), anEntry);
3010     if (argi > 1) aPartsDescr += ", ";
3011     aPartsDescr += anEntry;
3012   }
3013
3014   //Compute the Block value
3015   try {
3016 #if OCC_VERSION_LARGE > 0x06010000
3017     OCC_CATCH_SIGNALS;
3018 #endif
3019     // 1. Explode compound on solids
3020     TopTools_MapOfShape mapShape;
3021     Standard_Integer nbSolids = 0;
3022     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
3023     for (; exp.More(); exp.Next()) {
3024       if (mapShape.Add(exp.Current())) {
3025         nbSolids++;
3026       }
3027     }
3028
3029     mapShape.Clear();
3030     Standard_Integer ind = 1;
3031     TopTools_Array1OfShape aSolids (1, nbSolids);
3032     TColStd_Array1OfInteger aNbParts (1, nbSolids);
3033     for (exp.Init(aBlockOrComp, TopAbs_SOLID); exp.More(); exp.Next(), ind++) {
3034       if (mapShape.Add(exp.Current())) {
3035         TopoDS_Shape aSolid = exp.Current();
3036         aSolids(ind) = aSolid;
3037         aNbParts(ind) = 0;
3038
3039         // 2. Define quantity of parts, contained in each solid
3040         TopTools_IndexedMapOfShape aSubShapes;
3041         TopExp::MapShapes(aSolid, aSubShapes);
3042         for (argi = 1; argi <= aLen; argi++) {
3043           if (aSubShapes.Contains(anArgs(argi))) {
3044             aNbParts(ind)++;
3045           }
3046         }
3047       }
3048     }
3049
3050     // 3. Define solid, containing maximum quantity of parts
3051     Standard_Integer maxNb = 0, nbFound = 0;
3052     TopoDS_Shape aShape;
3053     for (ind = 1; ind <= nbSolids; ind++) {
3054       if (aNbParts(ind) > maxNb) {
3055         maxNb = aNbParts(ind);
3056         aShape = aSolids(ind);
3057         nbFound = 1;
3058       } else if (aNbParts(ind) == maxNb) {
3059         nbFound++;
3060       } else {
3061       }
3062     }
3063     if (nbFound > 1) {
3064       SetErrorCode("Multiple blocks, containing maximum quantity of the given parts, are found");
3065       return NULL;
3066     } else if (nbFound == 0) {
3067       SetErrorCode("There are no blocks, containing the given parts");
3068       return NULL;
3069     } else {
3070       TopTools_IndexedMapOfShape anIndices;
3071       TopExp::MapShapes(aBlockOrComp, anIndices);
3072       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
3073       anArray->SetValue(1, anIndices.FindIndex(aShape));
3074       aResult = GetEngine()->AddSubShape(theCompound, anArray);
3075     }
3076   } catch (Standard_Failure) {
3077     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3078     SetErrorCode(aFail->GetMessageString());
3079     return NULL;
3080   }
3081
3082   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
3083
3084   //Make a Python command
3085   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetBlockByParts("
3086     << theCompound << ", [" << aPartsDescr.ToCString() << "])";
3087
3088   SetErrorCode(OK);
3089   return aResult;
3090 }
3091
3092 //=============================================================================
3093 /*!
3094  *  GetBlocksByParts
3095  */
3096 //=============================================================================
3097 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::GetBlocksByParts
3098                       (Handle(GEOM_Object)                         theCompound,
3099                        const Handle(TColStd_HSequenceOfTransient)& theParts)
3100 {
3101   SetErrorCode(KO);
3102
3103   if (theCompound.IsNull() || theParts.IsNull()) return NULL;
3104   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
3105   if (aBlockOrComp.IsNull()) return NULL;
3106
3107   Handle(TColStd_HSequenceOfTransient) aBlocks = new TColStd_HSequenceOfTransient;
3108   Handle(GEOM_Object) anObj;
3109   Handle(GEOM_Function) aFunction;
3110
3111   //Get the parts
3112   Standard_Integer argi, aLen = theParts->Length();
3113   TopTools_Array1OfShape anArgs (1, aLen);
3114   TCollection_AsciiString anEntry, aPartsDescr, anAsciiList;
3115
3116   for (argi = 1; argi <= aLen; argi++) {
3117     Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast(theParts->Value(argi));
3118     Handle(GEOM_Function) aRef = anObj->GetLastFunction();
3119     if (aRef.IsNull()) return NULL;
3120
3121     TopoDS_Shape anArg = aRef->GetValue();
3122     if (anArg.IsNull()) {
3123       SetErrorCode("Null shape is given as argument");
3124       return NULL;
3125     }
3126     anArgs(argi) = anArg;
3127
3128     // For Python command
3129     TDF_Tool::Entry(anObj->GetEntry(), anEntry);
3130     aPartsDescr += anEntry + ", ";
3131   }
3132
3133   //Get the Blocks
3134   try {
3135 #if OCC_VERSION_LARGE > 0x06010000
3136     OCC_CATCH_SIGNALS;
3137 #endif
3138     TopTools_MapOfShape mapShape;
3139     Standard_Integer nbSolids = 0;
3140     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
3141     for (; exp.More(); exp.Next()) {
3142       if (mapShape.Add(exp.Current())) {
3143         nbSolids++;
3144       }
3145     }
3146
3147     mapShape.Clear();
3148     Standard_Integer ind = 1;
3149     TopTools_Array1OfShape aSolids (1, nbSolids);
3150     TColStd_Array1OfInteger aNbParts (1, nbSolids);
3151     for (exp.Init(aBlockOrComp, TopAbs_SOLID); exp.More(); exp.Next(), ind++) {
3152       if (mapShape.Add(exp.Current())) {
3153         TopoDS_Shape aSolid = exp.Current();
3154         aSolids(ind) = aSolid;
3155         aNbParts(ind) = 0;
3156
3157         // 2. Define quantity of parts, contained in each solid
3158         TopTools_IndexedMapOfShape aSubShapes;
3159         TopExp::MapShapes(aSolid, aSubShapes);
3160         for (argi = 1; argi <= aLen; argi++) {
3161           if (aSubShapes.Contains(anArgs(argi))) {
3162             aNbParts(ind)++;
3163           }
3164         }
3165       }
3166     }
3167
3168     // 3. Define solid, containing maximum quantity of parts
3169     Standard_Integer maxNb = 0, nbFound = 0;
3170     for (ind = 1; ind <= nbSolids; ind++) {
3171       if (aNbParts(ind) > maxNb) {
3172         maxNb = aNbParts(ind);
3173         nbFound = 1;
3174       } else if (aNbParts(ind) == maxNb) {
3175         nbFound++;
3176       } else {
3177       }
3178     }
3179     if (nbFound == 0) {
3180       SetErrorCode("There are no blocks, containing the given parts");
3181       return NULL;
3182     }
3183
3184     // Map shapes
3185     TopTools_IndexedMapOfShape anIndices;
3186     TopExp::MapShapes(aBlockOrComp, anIndices);
3187     Handle(TColStd_HArray1OfInteger) anArray;
3188
3189     for (ind = 1; ind <= nbSolids; ind++) {
3190       if (aNbParts(ind) == maxNb) {
3191         anArray = new TColStd_HArray1OfInteger(1,1);
3192         anArray->SetValue(1, anIndices.FindIndex(aSolids(ind)));
3193         anObj = GetEngine()->AddSubShape(theCompound, anArray);
3194         aBlocks->Append(anObj);
3195
3196         // For Python command
3197         TDF_Tool::Entry(anObj->GetEntry(), anEntry);
3198         anAsciiList += anEntry + ", ";
3199         if (aFunction.IsNull())
3200           aFunction = anObj->GetLastFunction();
3201       }
3202     }
3203   }
3204   catch (Standard_Failure) {
3205     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3206     SetErrorCode(aFail->GetMessageString());
3207     return NULL;
3208   }
3209
3210   //Make a Python command
3211   aPartsDescr.Trunc(aPartsDescr.Length() - 2);
3212   anAsciiList.Trunc(anAsciiList.Length() - 2);
3213
3214   GEOM::TPythonDump(aFunction) << "[" << anAsciiList.ToCString()
3215     << "] = geompy.GetBlocksByParts(" << theCompound
3216       << ", [" << aPartsDescr.ToCString() << "])";
3217
3218   SetErrorCode(OK);
3219   return aBlocks;
3220 }
3221
3222 //=============================================================================
3223 /*!
3224  *  MakeMultiTransformation1D
3225  */
3226 //=============================================================================
3227 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation1D
3228                                              (Handle(GEOM_Object)    theObject,
3229                                               const Standard_Integer theDirFace1,
3230                                               const Standard_Integer theDirFace2,
3231                                               const Standard_Integer theNbTimes)
3232 {
3233   SetErrorCode(KO);
3234
3235   if (theObject.IsNull()) return NULL;
3236
3237   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
3238   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be moved
3239
3240   //Add a new Copy object
3241   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
3242
3243   //Add a translate function
3244   Handle(GEOM_Function) aFunction =
3245     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_MULTI_TRANSFORM_1D);
3246
3247   //Check if the function is set correctly
3248   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
3249
3250   GEOMImpl_IBlockTrsf aTI (aFunction);
3251   aTI.SetOriginal(aLastFunction);
3252   aTI.SetFace1U(theDirFace1);
3253   aTI.SetFace2U(theDirFace2);
3254   aTI.SetNbIterU(theNbTimes);
3255
3256   //Compute the transformation
3257   try {
3258 #if OCC_VERSION_LARGE > 0x06010000
3259     OCC_CATCH_SIGNALS;
3260 #endif
3261     if (!GetSolver()->ComputeFunction(aFunction)) {
3262       SetErrorCode("Block driver failed to make multi-transformation");
3263       return NULL;
3264     }
3265   }
3266   catch (Standard_Failure) {
3267     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3268     SetErrorCode(aFail->GetMessageString());
3269     return NULL;
3270   }
3271
3272   //Make a Python command
3273   GEOM::TPythonDump(aFunction) << aCopy << " = geompy.MakeMultiTransformation1D("
3274     << theObject << ", " << theDirFace1 << ", " << theDirFace2 << ", " << theNbTimes << ")";
3275
3276   SetErrorCode(OK);
3277   return aCopy;
3278 }
3279
3280 //=============================================================================
3281 /*!
3282  *  MakeMultiTransformation2D
3283  */
3284 //=============================================================================
3285 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation2D
3286                                              (Handle(GEOM_Object)    theObject,
3287                                               const Standard_Integer theDirFace1U,
3288                                               const Standard_Integer theDirFace2U,
3289                                               const Standard_Integer theNbTimesU,
3290                                               const Standard_Integer theDirFace1V,
3291                                               const Standard_Integer theDirFace2V,
3292                                               const Standard_Integer theNbTimesV)
3293 {
3294   SetErrorCode(KO);
3295
3296   if (theObject.IsNull()) return NULL;
3297
3298   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
3299   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be moved
3300
3301   //Add a new Copy object
3302   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
3303
3304   //Add a translate function
3305   Handle(GEOM_Function) aFunction =
3306     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_MULTI_TRANSFORM_2D);
3307
3308   //Check if the function is set correctly
3309   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
3310
3311   GEOMImpl_IBlockTrsf aTI (aFunction);
3312   aTI.SetOriginal(aLastFunction);
3313   aTI.SetFace1U(theDirFace1U);
3314   aTI.SetFace2U(theDirFace2U);
3315   aTI.SetNbIterU(theNbTimesU);
3316   aTI.SetFace1V(theDirFace1V);
3317   aTI.SetFace2V(theDirFace2V);
3318   aTI.SetNbIterV(theNbTimesV);
3319
3320   //Compute the transformation
3321   try {
3322 #if OCC_VERSION_LARGE > 0x06010000
3323     OCC_CATCH_SIGNALS;
3324 #endif
3325     if (!GetSolver()->ComputeFunction(aFunction)) {
3326       SetErrorCode("Block driver failed to make multi-transformation");
3327       return NULL;
3328     }
3329   }
3330   catch (Standard_Failure) {
3331     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3332     SetErrorCode(aFail->GetMessageString());
3333     return NULL;
3334   }
3335
3336   //Make a Python command
3337   GEOM::TPythonDump(aFunction) << aCopy << " = geompy.MakeMultiTransformation2D("
3338     << theObject << ", " << theDirFace1U << ", " << theDirFace2U << ", " << theNbTimesU
3339       << ", " << theDirFace1V << ", " << theDirFace2V << ", " << theNbTimesV << ")";
3340
3341   SetErrorCode(OK);
3342   return aCopy;
3343 }
3344
3345 //=============================================================================
3346 /*!
3347  *  Propagate
3348  */
3349 //=============================================================================
3350 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::Propagate
3351                                                  (Handle(GEOM_Object) theShape)
3352 {
3353   SetErrorCode(KO);
3354
3355   if (theShape.IsNull()) return NULL;
3356
3357   TopoDS_Shape aShape = theShape->GetValue();
3358   if (aShape.IsNull()) return NULL;
3359
3360   TopTools_IndexedMapOfShape anIndices;
3361   TopExp::MapShapes(aShape, anIndices);
3362
3363   TopTools_IndexedDataMapOfShapeListOfShape MEW;
3364   GEOMImpl_Block6Explorer::MapShapesAndAncestors
3365     (aShape, TopAbs_EDGE, TopAbs_WIRE, MEW);
3366   Standard_Integer ie, nbEdges = MEW.Extent();
3367
3368   // Result
3369   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
3370
3371   TopTools_MapOfShape mapAcceptedEdges;
3372   TCollection_AsciiString aListRes, anEntry;
3373
3374   // Sort shapes in current chain (Mantis issue 21053)
3375   TopTools_DataMapOfShapeListOfShape aMapChains;
3376   TopTools_ListOfShape aFirstInChains;
3377
3378   for (ie = 1; ie <= nbEdges; ie++) {
3379     TopoDS_Shape curE = MEW.FindKey(ie);
3380
3381     if (mapAcceptedEdges.Contains(curE)) continue;
3382
3383     // Build the chain
3384     TopTools_ListOfShape currentChain;
3385     TopTools_ListOfShape listPrevEdges;
3386
3387     currentChain.Append(curE);
3388     listPrevEdges.Append(curE);
3389     mapAcceptedEdges.Add(curE);
3390
3391     // Collect all edges pass by pass
3392     while (listPrevEdges.Extent() > 0) {
3393       // List of edges, added to chain on this cycle pass
3394       TopTools_ListOfShape listCurEdges;
3395
3396       // Find the next portion of edges
3397       TopTools_ListIteratorOfListOfShape itE (listPrevEdges);
3398       for (; itE.More(); itE.Next()) {
3399         TopoDS_Shape anE = itE.Value();
3400
3401         // Iterate on faces, having edge <anE>
3402         TopTools_ListIteratorOfListOfShape itW (MEW.FindFromKey(anE));
3403         for (; itW.More(); itW.Next()) {
3404           TopoDS_Shape aW = itW.Value();
3405           TopoDS_Shape anOppE;
3406
3407           BRepTools_WireExplorer aWE (TopoDS::Wire(aW));
3408           Standard_Integer nb = 1, found = 0;
3409           TopTools_Array1OfShape anEdges (1,4);
3410           for (; aWE.More(); aWE.Next(), nb++) {
3411             if (nb > 4) {
3412               found = 0;
3413               break;
3414             }
3415             anEdges(nb) = aWE.Current();
3416             if (anEdges(nb).IsSame(anE)) found = nb;
3417           }
3418
3419           if (nb == 5 && found > 0) {
3420             // Quadrangle face found, get an opposite edge
3421             Standard_Integer opp = found + 2;
3422             if (opp > 4) opp -= 4;
3423             anOppE = anEdges(opp);
3424
3425             if (!mapAcceptedEdges.Contains(anOppE)) {
3426               // Add found edge to the chain
3427               currentChain.Append(anOppE);
3428               listCurEdges.Append(anOppE);
3429               mapAcceptedEdges.Add(anOppE);
3430             }
3431           } // if (nb == 5 && found > 0)
3432         } // for (; itF.More(); itF.Next())
3433       } // for (; itE.More(); itE.Next())
3434
3435       listPrevEdges = listCurEdges;
3436     } // while (listPrevEdges.Extent() > 0)
3437
3438     // Sort shapes in current chain (Mantis issue 21053)
3439     GEOMUtils::SortShapes(currentChain, Standard_False);
3440     aFirstInChains.Append(currentChain.First());
3441     aMapChains.Bind(currentChain.First(), currentChain);
3442   }
3443
3444   // Sort chains (Mantis issue 21053)
3445   GEOMUtils::SortShapes(aFirstInChains, Standard_False);
3446
3447   // Store sorted chains in the document
3448   TopTools_ListIteratorOfListOfShape aChainsIt (aFirstInChains);
3449   for (; aChainsIt.More(); aChainsIt.Next()) {
3450     TopoDS_Shape aFirstInChain = aChainsIt.Value();
3451     const TopTools_ListOfShape& currentChain = aMapChains.Find(aFirstInChain);
3452
3453     // Store the chain in the document
3454     Handle(TColStd_HArray1OfInteger) anArray =
3455       new TColStd_HArray1OfInteger (1, currentChain.Extent());
3456
3457     // Fill array of sub-shape indices
3458     TopTools_ListIteratorOfListOfShape itSub (currentChain);
3459     for (int index = 1; itSub.More(); itSub.Next(), ++index) {
3460       int id = anIndices.FindIndex(itSub.Value());
3461       anArray->SetValue(index, id);
3462     }
3463
3464     // Add a new group object
3465     Handle(GEOM_Object) aChain = GetEngine()->AddSubShape(theShape, anArray);
3466
3467     // Set a GROUP type
3468     aChain->SetType(GEOM_GROUP);
3469
3470     // Set a sub-shape type
3471     TDF_Label aFreeLabel = aChain->GetFreeLabel();
3472     TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)TopAbs_EDGE);
3473
3474     // Add the chain to the result
3475     aSeq->Append(aChain);
3476
3477     //Make a Python command
3478     TDF_Tool::Entry(aChain->GetEntry(), anEntry);
3479     aListRes += anEntry + ", ";
3480   }
3481
3482   if (aSeq->IsEmpty()) {
3483     SetErrorCode("There are no quadrangle faces in the shape");
3484     return aSeq;
3485   }
3486
3487   aListRes.Trunc(aListRes.Length() - 2);
3488
3489   // The Propagation doesn't change object so no new function is required.
3490   Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
3491
3492   // Make a Python command
3493   GEOM::TPythonDump(aFunction, /*append=*/true)
3494     << "[" << aListRes.ToCString() << "] = geompy.Propagate(" << theShape << ")";
3495
3496   SetErrorCode(OK);
3497   return aSeq;
3498 }