src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx

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