src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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