src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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