src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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