src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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