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