src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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