src/GEOMImpl/GEOMImpl_IShapesOperations.cxx

   1 // Copyright (C) 2005  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   2 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   3 //
   4 // This library is free software; you can redistribute it and/or
   5 // modify it under the terms of the GNU Lesser General Public
   6 // License as published by the Free Software Foundation; either
   7 // version 2.1 of the License.
   8 //
   9 // This library is distributed in the hope that it will be useful
  10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 // Lesser General Public License for more details.
  13 //
  14 // You should have received a copy of the GNU Lesser General Public
  15 // License along with this library; if not, write to the Free Software
  16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 //
  18 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 //
  20 // File      : GEOMImpl_IShapesOperations.cxx
  21 // Created   :
  22 // Author    : modified by Lioka RAZAFINDRAZAKA (CEA) 22/06/2007
  23 // Project   : SALOME
  24 // $Header$
  25
  26 #include <Standard_Stream.hxx>
  27
  28 #include "GEOMImpl_IShapesOperations.hxx"
  29
  30 #include "GEOMImpl_Types.hxx"
  31
  32 #include "GEOMImpl_VectorDriver.hxx"
  33 #include "GEOMImpl_ShapeDriver.hxx"
  34 #include "GEOMImpl_CopyDriver.hxx"
  35 #include "GEOMImpl_GlueDriver.hxx"
  36
  37 #include "GEOMImpl_IVector.hxx"
  38 #include "GEOMImpl_IShapes.hxx"
  39 #include "GEOMImpl_IGlue.hxx"
  40
  41 #include "GEOMImpl_Block6Explorer.hxx"
  42
  43 #include "GEOM_Function.hxx"
  44 #include "GEOM_ISubShape.hxx"
  45 #include "GEOM_PythonDump.hxx"
  46
  47 #include "GEOMAlgo_FinderShapeOn1.hxx"
  48 #include "GEOMAlgo_FinderShapeOnQuad.hxx"
  49 #include "GEOMAlgo_FinderShapeOn2.hxx"
  50 #include "GEOMAlgo_ClsfBox.hxx"
  51 #include "GEOMAlgo_ClsfSolid.hxx"
  52 #include "GEOMAlgo_Gluer1.hxx"
  53 #include "GEOMAlgo_ListIteratorOfListOfCoupleOfShapes.hxx"
  54 #include "GEOMAlgo_CoupleOfShapes.hxx"
  55 #include "GEOMAlgo_ListOfCoupleOfShapes.hxx"
  56
  57 #include "utilities.h"
  58 #include "OpUtil.hxx"
  59 #include "Utils_ExceptHandlers.hxx"
  60
  61 #include <TFunction_DriverTable.hxx>
  62 #include <TFunction_Driver.hxx>
  63 #include <TFunction_Logbook.hxx>
  64 #include <TDataStd_Integer.hxx>
  65 #include <TDataStd_IntegerArray.hxx>
  66 #include <TDF_Tool.hxx>
  67
  68 #include <BRepExtrema_ExtCF.hxx>
  69 #include <BRepExtrema_DistShapeShape.hxx>
  70
  71 #include <BRep_Tool.hxx>
  72 #include <BRep_Builder.hxx>
  73 #include <BRepTools.hxx>
  74 #include <BRepGProp.hxx>
  75 #include <BRepAdaptor_Curve.hxx>
  76 #include <BRepAdaptor_Surface.hxx>
  77 #include <BRepBndLib.hxx>
  78 #include <BRepBuilderAPI_MakeFace.hxx>
  79 #include <BRepMesh_IncrementalMesh.hxx>
  80
  81 #include <TopAbs.hxx>
  82 #include <TopExp.hxx>
  83 #include <TopoDS.hxx>
  84 #include <TopoDS_Shape.hxx>
  85 #include <TopoDS_Solid.hxx>
  86 #include <TopoDS_Face.hxx>
  87 #include <TopoDS_Edge.hxx>
  88 #include <TopoDS_Vertex.hxx>
  89 #include <TopoDS_Compound.hxx>
  90 #include <TopoDS_Iterator.hxx>
  91 #include <TopExp_Explorer.hxx>
  92 #include <TopLoc_Location.hxx>
  93 #include <TopTools_MapOfShape.hxx>
  94 #include <TopTools_MapOfOrientedShape.hxx>
  95 #include <TopTools_Array1OfShape.hxx>
  96 #include <TopTools_ListIteratorOfListOfShape.hxx>
  97 #include <TopTools_IndexedMapOfShape.hxx>
  98
  99 #include <Geom_Surface.hxx>
 100 #include <Geom_Plane.hxx>
 101 #include <Geom_SphericalSurface.hxx>
 102 #include <Geom_CylindricalSurface.hxx>
 103 #include <GeomAdaptor_Surface.hxx>
 104
 105 #include <GeomLib_Tool.hxx>
 106 #include <Geom2d_Curve.hxx>
 107
 108 #include <Bnd_Box.hxx>
 109 #include <GProp_GProps.hxx>
 110 #include <gp_Pnt.hxx>
 111 #include <gp_Lin.hxx>
 112 #include <TColStd_ListOfInteger.hxx>
 113 #include <TColStd_ListIteratorOfListOfInteger.hxx>
 114 #include <TColStd_Array1OfReal.hxx>
 115 #include <TColStd_HArray1OfInteger.hxx>
 116
 117 #include <vector>
 118
 119 #include <Standard_NullObject.hxx>
 120 #include <Standard_Failure.hxx>
 121 #include <Standard_ErrorHandler.hxx> // CAREFUL ! position of this file is critic : see Lucien PIGNOLONI / OCC
 122
 123 // Includes added for GetInPlace algorithm improvement
 124
 125 #include <GEOMImpl_MeasureDriver.hxx>
 126 #include <GEOMImpl_IMeasure.hxx>
 127 #include <BRepBuilderAPI_MakeVertex.hxx>
 128
 129 #include <BRepClass_FaceClassifier.hxx>
 130 #include <BRepClass3d_SolidClassifier.hxx>
 131 #include <Precision.hxx>
 132
 133 //=============================================================================
 134 /*!
 135  *   constructor:
 136  */
 137 //=============================================================================
 138 GEOMImpl_IShapesOperations::GEOMImpl_IShapesOperations (GEOM_Engine* theEngine, int theDocID)
 139 : GEOM_IOperations(theEngine, theDocID)
 140 {
 141   MESSAGE("GEOMImpl_IShapesOperations::GEOMImpl_IShapesOperations");
 142 }
 143
 144 //=============================================================================
 145 /*!
 146  *  destructor
 147  */
 148 //=============================================================================
 149 GEOMImpl_IShapesOperations::~GEOMImpl_IShapesOperations()
 150 {
 151   MESSAGE("GEOMImpl_IShapesOperations::~GEOMImpl_IShapesOperations");
 152 }
 153
 154
 155 //=============================================================================
 156 /*!
 157  *  MakeEdge
 158  */
 159 //=============================================================================
 160 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeEdge
 161                      (Handle(GEOM_Object) thePnt1, Handle(GEOM_Object) thePnt2)
 162 {
 163   SetErrorCode(KO);
 164
 165   if (thePnt1.IsNull() || thePnt2.IsNull()) return NULL;
 166
 167   //Add a new Edge object
 168   Handle(GEOM_Object) anEdge = GetEngine()->AddObject(GetDocID(), GEOM_EDGE);
 169
 170   //Add a new Vector function
 171   Handle(GEOM_Function) aFunction =
 172     anEdge->AddFunction(GEOMImpl_VectorDriver::GetID(), VECTOR_TWO_PNT);
 173
 174   //Check if the function is set correctly
 175   if (aFunction->GetDriverGUID() != GEOMImpl_VectorDriver::GetID()) return NULL;
 176
 177   GEOMImpl_IVector aPI (aFunction);
 178
 179   Handle(GEOM_Function) aRef1 = thePnt1->GetLastFunction();
 180   Handle(GEOM_Function) aRef2 = thePnt2->GetLastFunction();
 181   if (aRef1.IsNull() || aRef2.IsNull()) return NULL;
 182
 183   aPI.SetPoint1(aRef1);
 184   aPI.SetPoint2(aRef2);
 185
 186   //Compute the Edge value
 187   try {
 188 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 189     OCC_CATCH_SIGNALS;
 190 #endif
 191     if (!GetSolver()->ComputeFunction(aFunction)) {
 192       SetErrorCode("Vector driver failed");
 193       return NULL;
 194     }
 195   }
 196   catch (Standard_Failure) {
 197     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 198     SetErrorCode(aFail->GetMessageString());
 199     return NULL;
 200   }
 201
 202   //Make a Python command
 203   GEOM::TPythonDump(aFunction) << anEdge << " = geompy.MakeEdge("
 204                                << thePnt1 << ", " << thePnt2 << ")";
 205
 206   SetErrorCode(OK);
 207   return anEdge;
 208 }
 209
 210 //=============================================================================
 211 /*!
 212  *  MakeWire
 213  */
 214 //=============================================================================
 215 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeWire
 216                              (std::list<Handle(GEOM_Object)> theShapes)
 217 {
 218   return MakeShape(theShapes, GEOM_WIRE, WIRE_EDGES, "MakeWire");
 219 }
 220
 221 //=============================================================================
 222 /*!
 223  *  MakeFace
 224  */
 225 //=============================================================================
 226 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeFace (Handle(GEOM_Object) theWire,
 227                                                           const bool isPlanarWanted)
 228 {
 229   SetErrorCode(KO);
 230
 231   if (theWire.IsNull()) return NULL;
 232
 233   //Add a new Face object
 234   Handle(GEOM_Object) aFace = GetEngine()->AddObject(GetDocID(), GEOM_FACE);
 235
 236   //Add a new Shape function for creation of a face from a wire
 237   Handle(GEOM_Function) aFunction =
 238     aFace->AddFunction(GEOMImpl_ShapeDriver::GetID(), FACE_WIRE);
 239   if (aFunction.IsNull()) return NULL;
 240
 241   //Check if the function is set correctly
 242   if (aFunction->GetDriverGUID() != GEOMImpl_ShapeDriver::GetID()) return NULL;
 243
 244   GEOMImpl_IShapes aCI (aFunction);
 245
 246   Handle(GEOM_Function) aRefWire = theWire->GetLastFunction();
 247
 248   if (aRefWire.IsNull()) return NULL;
 249
 250   aCI.SetBase(aRefWire);
 251   aCI.SetIsPlanar(isPlanarWanted);
 252
 253   //Compute the Face value
 254   try {
 255 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 256     OCC_CATCH_SIGNALS;
 257 #endif
 258     if (!GetSolver()->ComputeFunction(aFunction)) {
 259       SetErrorCode("Shape driver failed to compute a face");
 260       return NULL;
 261     }
 262   }
 263   catch (Standard_Failure) {
 264     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 265     SetErrorCode(aFail->GetMessageString());
 266     return NULL;
 267   }
 268
 269   //Make a Python command
 270   GEOM::TPythonDump(aFunction) << aFace << " = geompy.MakeFace("
 271     << theWire << ", " << (int)isPlanarWanted << ")";
 272
 273   SetErrorCode(OK);
 274   return aFace;
 275 }
 276
 277 //=============================================================================
 278 /*!
 279  *  MakeFaceWires
 280  */
 281 //=============================================================================
 282 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeFaceWires
 283                              (std::list<Handle(GEOM_Object)> theShapes,
 284                               const bool isPlanarWanted)
 285 {
 286   SetErrorCode(KO);
 287
 288   //Add a new object
 289   Handle(GEOM_Object) aShape = GetEngine()->AddObject(GetDocID(), GEOM_FACE);
 290
 291   //Add a new function
 292   Handle(GEOM_Function) aFunction =
 293     aShape->AddFunction(GEOMImpl_ShapeDriver::GetID(), FACE_WIRES);
 294   if (aFunction.IsNull()) return NULL;
 295
 296   //Check if the function is set correctly
 297   if (aFunction->GetDriverGUID() != GEOMImpl_ShapeDriver::GetID()) return NULL;
 298
 299   GEOMImpl_IShapes aCI (aFunction);
 300
 301   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
 302
 303   // Shapes
 304   std::list<Handle(GEOM_Object)>::iterator it = theShapes.begin();
 305   for (; it != theShapes.end(); it++) {
 306     Handle(GEOM_Function) aRefSh = (*it)->GetLastFunction();
 307     if (aRefSh.IsNull()) {
 308       SetErrorCode("NULL argument shape for the face construction");
 309       return NULL;
 310     }
 311     aShapesSeq->Append(aRefSh);
 312   }
 313   aCI.SetShapes(aShapesSeq);
 314
 315   aCI.SetIsPlanar(isPlanarWanted);
 316
 317   //Compute the shape
 318   try {
 319 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 320     OCC_CATCH_SIGNALS;
 321 #endif
 322     if (!GetSolver()->ComputeFunction(aFunction)) {
 323       SetErrorCode("Shape driver failed");
 324       return NULL;
 325     }
 326   }
 327   catch (Standard_Failure) {
 328     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 329     SetErrorCode(aFail->GetMessageString());
 330     return NULL;
 331   }
 332
 333   //Make a Python command
 334   GEOM::TPythonDump pd (aFunction);
 335   pd << aShape << " = geompy.MakeFaceWires([";
 336
 337   // Shapes
 338   it = theShapes.begin();
 339   if (it != theShapes.end()) {
 340     pd << (*it++);
 341     while (it != theShapes.end()) {
 342       pd << ", " << (*it++);
 343     }
 344   }
 345   pd << "], " << (int)isPlanarWanted << ")";
 346
 347   SetErrorCode(OK);
 348   return aShape;
 349 }
 350
 351 //=============================================================================
 352 /*!
 353  *  MakeShell
 354  */
 355 //=============================================================================
 356 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeShell
 357                              (std::list<Handle(GEOM_Object)> theShapes)
 358 {
 359   return MakeShape(theShapes, GEOM_SHELL, SHELL_FACES, "MakeShell");
 360 }
 361
 362 //=============================================================================
 363 /*!
 364  *  MakeSolidShells
 365  */
 366 //=============================================================================
 367 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeSolidShells
 368                              (std::list<Handle(GEOM_Object)> theShapes)
 369 {
 370   return MakeShape(theShapes, GEOM_SOLID, SOLID_SHELLS, "MakeSolid");
 371 }
 372
 373 //=============================================================================
 374 /*!
 375  *  MakeSolidShell
 376  */
 377 //=============================================================================
 378 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeSolidShell (Handle(GEOM_Object) theShell)
 379 {
 380   SetErrorCode(KO);
 381
 382   if (theShell.IsNull()) return NULL;
 383
 384   //Add a new Solid object
 385   Handle(GEOM_Object) aSolid = GetEngine()->AddObject(GetDocID(), GEOM_SOLID);
 386
 387   //Add a new Solid function for creation of a solid from a shell
 388   Handle(GEOM_Function) aFunction =
 389     aSolid->AddFunction(GEOMImpl_ShapeDriver::GetID(), SOLID_SHELL);
 390   if (aFunction.IsNull()) return NULL;
 391
 392   //Check if the function is set correctly
 393   if (aFunction->GetDriverGUID() != GEOMImpl_ShapeDriver::GetID()) return NULL;
 394
 395   GEOMImpl_IShapes aCI (aFunction);
 396
 397   Handle(GEOM_Function) aRefShell = theShell->GetLastFunction();
 398
 399   if (aRefShell.IsNull()) return NULL;
 400
 401   aCI.SetBase(aRefShell);
 402
 403   //Compute the Solid value
 404   try {
 405 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 406     OCC_CATCH_SIGNALS;
 407 #endif
 408     if (!GetSolver()->ComputeFunction(aFunction)) {
 409       SetErrorCode("Solid driver failed");
 410       return NULL;
 411     }
 412   }
 413   catch (Standard_Failure) {
 414     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 415     SetErrorCode(aFail->GetMessageString());
 416     return NULL;
 417   }
 418
 419   //Make a Python command
 420   GEOM::TPythonDump(aFunction) << aSolid
 421     << " = geompy.MakeSolid(" << theShell << ")";
 422
 423   SetErrorCode(OK);
 424   return aSolid;
 425 }
 426
 427 //=============================================================================
 428 /*!
 429  *  MakeCompound
 430  */
 431 //=============================================================================
 432 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeCompound
 433                              (std::list<Handle(GEOM_Object)> theShapes)
 434 {
 435   return MakeShape(theShapes, GEOM_COMPOUND, COMPOUND_SHAPES, "MakeCompound");
 436 }
 437
 438 //=============================================================================
 439 /*!
 440  *  MakeShape
 441  */
 442 //=============================================================================
 443 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeShape
 444                              (std::list<Handle(GEOM_Object)>      theShapes,
 445                               const Standard_Integer         theObjectType,
 446                               const Standard_Integer         theFunctionType,
 447                               const TCollection_AsciiString& theMethodName)
 448 {
 449   SetErrorCode(KO);
 450
 451   //Add a new object
 452   Handle(GEOM_Object) aShape = GetEngine()->AddObject(GetDocID(), theObjectType);
 453
 454   //Add a new function
 455   Handle(GEOM_Function) aFunction =
 456     aShape->AddFunction(GEOMImpl_ShapeDriver::GetID(), theFunctionType);
 457   if (aFunction.IsNull()) return NULL;
 458
 459   //Check if the function is set correctly
 460   if (aFunction->GetDriverGUID() != GEOMImpl_ShapeDriver::GetID()) return NULL;
 461
 462   GEOMImpl_IShapes aCI (aFunction);
 463
 464   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
 465
 466   // Shapes
 467   std::list<Handle(GEOM_Object)>::iterator it = theShapes.begin();
 468   for (; it != theShapes.end(); it++) {
 469     Handle(GEOM_Function) aRefSh = (*it)->GetLastFunction();
 470     if (aRefSh.IsNull()) {
 471       SetErrorCode("NULL argument shape for the shape construction");
 472       return NULL;
 473     }
 474     aShapesSeq->Append(aRefSh);
 475   }
 476   aCI.SetShapes(aShapesSeq);
 477
 478   //Compute the shape
 479   try {
 480 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 481     OCC_CATCH_SIGNALS;
 482 #endif
 483     if (!GetSolver()->ComputeFunction(aFunction)) {
 484       SetErrorCode("Shape driver failed");
 485       return NULL;
 486     }
 487   }
 488   catch (Standard_Failure) {
 489     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 490     SetErrorCode(aFail->GetMessageString());
 491     return NULL;
 492   }
 493
 494   //Make a Python command
 495   GEOM::TPythonDump pd (aFunction);
 496   pd << aShape << " = geompy." << theMethodName.ToCString() << "([";
 497
 498   // Shapes
 499   it = theShapes.begin();
 500   if (it != theShapes.end()) {
 501     pd << (*it++);
 502     while (it != theShapes.end()) {
 503       pd << ", " << (*it++);
 504     }
 505   }
 506   pd << "])";
 507
 508   SetErrorCode(OK);
 509   return aShape;
 510 }
 511
 512 //=============================================================================
 513 /*!
 514  *  MakeGlueFaces
 515  */
 516 //=============================================================================
 517 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeGlueFaces
 518                                                 (Handle(GEOM_Object) theShape,
 519                                                  const Standard_Real theTolerance,
 520                                                  const Standard_Boolean doKeepNonSolids)
 521 {
 522   SetErrorCode(KO);
 523
 524   if (theShape.IsNull()) return NULL;
 525
 526   //Add a new Glued object
 527   Handle(GEOM_Object) aGlued = GetEngine()->AddObject(GetDocID(), GEOM_GLUED);
 528
 529   //Add a new Glue function
 530   Handle(GEOM_Function) aFunction;
 531   aFunction = aGlued->AddFunction(GEOMImpl_GlueDriver::GetID(), GLUE_FACES);
 532   if (aFunction.IsNull()) return NULL;
 533
 534   //Check if the function is set correctly
 535   if (aFunction->GetDriverGUID() != GEOMImpl_GlueDriver::GetID()) return NULL;
 536
 537   GEOMImpl_IGlue aCI (aFunction);
 538
 539   Handle(GEOM_Function) aRefShape = theShape->GetLastFunction();
 540   if (aRefShape.IsNull()) return NULL;
 541
 542   aCI.SetBase(aRefShape);
 543   aCI.SetTolerance(theTolerance);
 544   aCI.SetKeepNonSolids(doKeepNonSolids);
 545
 546   //Compute the sub-shape value
 547   Standard_Boolean isWarning = Standard_False;
 548   try {
 549 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 550     OCC_CATCH_SIGNALS;
 551 #endif
 552     if (!GetSolver()->ComputeFunction(aFunction)) {
 553       SetErrorCode("Shape driver failed to glue faces");
 554       return NULL;
 555     }
 556   }
 557   catch (Standard_Failure) {
 558     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 559     SetErrorCode(aFail->GetMessageString());
 560     // to provide warning
 561     if (!aFunction->GetValue().IsNull()) {
 562       isWarning = Standard_True;
 563     } else {
 564       return NULL;
 565     }
 566   }
 567
 568   //Make a Python command
 569   GEOM::TPythonDump(aFunction) << aGlued << " = geompy.MakeGlueFaces("
 570     << theShape << ", " << theTolerance << ")";
 571
 572   // to provide warning
 573   if (!isWarning) SetErrorCode(OK);
 574   return aGlued;
 575 }
 576
 577 //=============================================================================
 578 /*!
 579  *  GetGlueFaces
 580  */
 581 //=============================================================================
 582 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetGlueFaces
 583                                                 (Handle(GEOM_Object) theShape,
 584                                                  const Standard_Real theTolerance)
 585 {
 586   SetErrorCode(KO);
 587
 588   if (theShape.IsNull()) return NULL;
 589   TopoDS_Shape aShape = theShape->GetValue();
 590   if (aShape.IsNull()) return NULL;
 591
 592   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
 593
 594   Standard_Integer iErr;
 595   TopoDS_Shape aS;
 596   GEOMAlgo_Gluer1 aGluer;
 597   GEOMAlgo_ListIteratorOfListOfCoupleOfShapes aItCS;
 598   GEOMAlgo_CoupleOfShapes aCS;
 599   GEOMAlgo_ListOfCoupleOfShapes aLCS;
 600
 601   //aGluer = new GEOMAlgo_Gluer1;
 602   aGluer.SetShape(aShape);
 603   aGluer.SetTolerance(theTolerance);
 604   aGluer.Perform();
 605   iErr = aGluer.ErrorStatus();
 606   if (iErr) return NULL;
 607
 608   TopTools_ListOfShape listShape;
 609   const GEOMAlgo_ListOfCoupleOfShapes& aLCSG = aGluer.GluedFaces();
 610   // Access to faces
 611   aItCS.Initialize(aLCSG);
 612   for (; aItCS.More(); aItCS.Next()) {
 613     const GEOMAlgo_CoupleOfShapes& aCSG = aItCS.Value();
 614     listShape.Append(aCSG.Shape1());
 615   }
 616
 617   TopTools_ListIteratorOfListOfShape itSub (listShape);
 618   TCollection_AsciiString anAsciiList, anEntry;
 619   TopTools_IndexedMapOfShape anIndices;
 620   TopExp::MapShapes(aShape, anIndices);
 621   Handle(TColStd_HArray1OfInteger) anArray;
 622   Handle(GEOM_Object) anObj;
 623   for (int index = 1; itSub.More(); itSub.Next(), ++index) {
 624     TopoDS_Shape aValue = itSub.Value();
 625     anArray = new TColStd_HArray1OfInteger(1,1);
 626     anArray->SetValue(1, anIndices.FindIndex(aValue));
 627     anObj = GetEngine()->AddSubShape(theShape, anArray);
 628     if (!anObj.IsNull()) {
 629       aSeq->Append(anObj);
 630
 631       // for python command
 632       TDF_Tool::Entry(anObj->GetEntry(), anEntry);
 633       anAsciiList += anEntry;
 634       anAsciiList += ",";
 635     }
 636   }
 637
 638   //Make a Python command
 639   if(anAsciiList.Length()>0)
 640     anAsciiList.Trunc(anAsciiList.Length() - 1);
 641   Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
 642   GEOM::TPythonDump pd (aFunction, /*append=*/true);
 643   pd << "[" << anAsciiList.ToCString();
 644   pd << "] = geompy.GetGlueFaces(" << theShape << ", " << theTolerance << ")";
 645
 646   SetErrorCode(OK);
 647
 648   return aSeq;
 649 }
 650
 651
 652 //=============================================================================
 653 /*!
 654  *  MakeGlueFacesByList
 655  */
 656 //=============================================================================
 657 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeGlueFacesByList
 658                                                 (Handle(GEOM_Object) theShape,
 659                                                  const Standard_Real theTolerance,
 660                                                  std::list<Handle(GEOM_Object)> theFaces,
 661                                                  const Standard_Boolean doKeepNonSolids)
 662 {
 663   SetErrorCode(KO);
 664
 665   if (theShape.IsNull()) return NULL;
 666
 667   //Add a new Glued object
 668   Handle(GEOM_Object) aGlued = GetEngine()->AddObject(GetDocID(), GEOM_GLUED);
 669
 670   //Add a new Glue function
 671   Handle(GEOM_Function) aFunction;
 672   aFunction = aGlued->AddFunction(GEOMImpl_GlueDriver::GetID(), GLUE_FACES_BY_LIST);
 673   if (aFunction.IsNull()) return NULL;
 674
 675   //Check if the function is set correctly
 676   if (aFunction->GetDriverGUID() != GEOMImpl_GlueDriver::GetID()) return NULL;
 677
 678   GEOMImpl_IGlue aCI (aFunction);
 679
 680   Handle(GEOM_Function) aRefShape = theShape->GetLastFunction();
 681   if (aRefShape.IsNull()) return NULL;
 682
 683   aCI.SetBase(aRefShape);
 684   aCI.SetTolerance(theTolerance);
 685   aCI.SetKeepNonSolids(doKeepNonSolids);
 686
 687   Handle(TColStd_HSequenceOfTransient) aFaces = new TColStd_HSequenceOfTransient;
 688   std::list<Handle(GEOM_Object)>::iterator it = theFaces.begin();
 689   for (; it != theFaces.end(); it++) {
 690     Handle(GEOM_Function) aRefSh = (*it)->GetLastFunction();
 691     if (aRefSh.IsNull()) {
 692       SetErrorCode("NULL argument shape for the shape construction");
 693       return NULL;
 694     }
 695     aFaces->Append(aRefSh);
 696   }
 697   aCI.SetFaces(aFaces);
 698
 699   //Compute the sub-shape value
 700   Standard_Boolean isWarning = Standard_False;
 701   try {
 702 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 703     OCC_CATCH_SIGNALS;
 704 #endif
 705     if (!GetSolver()->ComputeFunction(aFunction)) {
 706       SetErrorCode("Shape driver failed to glue faces");
 707       return NULL;
 708     }
 709   }
 710   catch (Standard_Failure) {
 711     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
 712     SetErrorCode(aFail->GetMessageString());
 713     // to provide warning
 714     if (!aFunction->GetValue().IsNull()) {
 715       isWarning = Standard_True;
 716     } else {
 717       return NULL;
 718     }
 719   }
 720
 721   //Make a Python command
 722
 723   GEOM::TPythonDump pd(aFunction);
 724   pd << aGlued << " = geompy.MakeGlueFacesByList("
 725      << theShape << ", " << theTolerance << ", [";
 726   // Faces
 727   it = theFaces.begin();
 728   if (it != theFaces.end()) {
 729     pd << (*it++);
 730     while (it != theFaces.end()) {
 731       pd << ", " << (*it++);
 732     }
 733   }
 734   pd << "])";
 735
 736
 737   // to provide warning
 738   if (!isWarning) SetErrorCode(OK);
 739   return aGlued;
 740 }
 741
 742
 743
 744 //=============================================================================
 745 /*!
 746  *  MakeExplode
 747  */
 748 //=============================================================================
 749 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::MakeExplode
 750                                           (Handle(GEOM_Object)    theShape,
 751                                            const Standard_Integer theShapeType,
 752                                            const Standard_Boolean isSorted)
 753 {
 754   SetErrorCode(KO);
 755
 756   if (theShape.IsNull()) return NULL;
 757   TopoDS_Shape aShape = theShape->GetValue();
 758   if (aShape.IsNull()) return NULL;
 759
 760   Handle(GEOM_Function) aMainShape = theShape->GetLastFunction();
 761
 762   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
 763   Handle(GEOM_Object) anObj;
 764   TopTools_MapOfShape mapShape;
 765   TopTools_ListOfShape listShape;
 766
 767   if (aShape.ShapeType() == TopAbs_COMPOUND &&
 768       (TopAbs_ShapeEnum(theShapeType) == TopAbs_SHAPE ||
 769        TopAbs_ShapeEnum(theShapeType) == TopAbs_COMPSOLID ||
 770        TopAbs_ShapeEnum(theShapeType) == TopAbs_COMPOUND)) {
 771     TopoDS_Iterator It (aShape, Standard_True, Standard_True);
 772     for (; It.More(); It.Next()) {
 773       if (mapShape.Add(It.Value())) {
 774         if (TopAbs_ShapeEnum(theShapeType) == TopAbs_SHAPE ||
 775             TopAbs_ShapeEnum(theShapeType) == It.Value().ShapeType()) {
 776           listShape.Append(It.Value());
 777         }
 778       }
 779     }
 780   } else {
 781     TopExp_Explorer exp (aShape, TopAbs_ShapeEnum(theShapeType));
 782     for (; exp.More(); exp.Next())
 783       if (mapShape.Add(exp.Current()))
 784         listShape.Append(exp.Current());
 785   }
 786
 787   if (listShape.IsEmpty()) {
 788     //SetErrorCode("The given shape has no sub-shapes of the requested type");
 789     SetErrorCode(NOT_FOUND_ANY); // NPAL18017
 790     return aSeq;
 791   }
 792
 793   if (isSorted)
 794     SortShapes(listShape);
 795
 796   TopTools_IndexedMapOfShape anIndices;
 797   TopExp::MapShapes(aShape, anIndices);
 798   Handle(TColStd_HArray1OfInteger) anArray;
 799
 800   TopTools_ListIteratorOfListOfShape itSub (listShape);
 801   TCollection_AsciiString anAsciiList, anEntry;
 802   for (int index = 1; itSub.More(); itSub.Next(), ++index)
 803   {
 804     TopoDS_Shape aValue = itSub.Value();
 805     anArray = new TColStd_HArray1OfInteger(1,1);
 806     anArray->SetValue(1, anIndices.FindIndex(aValue));
 807
 808     //anObj = GetEngine()->AddSubShape(theShape, anArray);
 809     {
 810       anObj = GetEngine()->AddObject(GetDocID(), GEOM_SUBSHAPE);
 811       Handle(GEOM_Function) aFunction = anObj->AddFunction(GEOM_Object::GetSubShapeID(), 1);
 812       if (aFunction.IsNull()) return aSeq;
 813
 814       GEOM_ISubShape aSSI (aFunction);
 815       aSSI.SetMainShape(aMainShape);
 816       aSSI.SetIndices(anArray);
 817
 818       // Set function value directly, as we know it.
 819       // Usage of Solver here would lead to significant loss of time,
 820       // because GEOM_SubShapeDriver will build TopTools_IndexedMapOfShape
 821       // on the main shape for each being calculated sub-shape separately.
 822       aFunction->SetValue(aValue);
 823     }
 824
 825     if (!anObj.IsNull()) {
 826       aSeq->Append(anObj);
 827
 828       // for python command
 829       TDF_Tool::Entry(anObj->GetEntry(), anEntry);
 830       anAsciiList += anEntry;
 831       anAsciiList += ",";
 832     }
 833   }
 834
 835   //Make a Python command
 836   anAsciiList.Trunc(anAsciiList.Length() - 1);
 837
 838   GEOM::TPythonDump pd (aMainShape, /*append=*/true);
 839   pd << "[" << anAsciiList.ToCString();
 840   pd << "] = geompy.SubShapeAll" << (isSorted ? "Sorted(" : "(");
 841   pd << theShape << ", " << TopAbs_ShapeEnum(theShapeType) << ")";
 842
 843   SetErrorCode(OK);
 844
 845   return aSeq;
 846 }
 847
 848 //=============================================================================
 849 /*!
 850  *  SubShapeAllIDs
 851  */
 852 //=============================================================================
 853 Handle(TColStd_HSequenceOfInteger) GEOMImpl_IShapesOperations::SubShapeAllIDs
 854                                           (Handle(GEOM_Object)    theShape,
 855                                            const Standard_Integer theShapeType,
 856                                            const Standard_Boolean isSorted)
 857 {
 858   SetErrorCode(KO);
 859
 860   if (theShape.IsNull()) return NULL;
 861   TopoDS_Shape aShape = theShape->GetValue();
 862   if (aShape.IsNull()) return NULL;
 863
 864   Handle(TColStd_HSequenceOfInteger) aSeq = new TColStd_HSequenceOfInteger;
 865   TopTools_MapOfShape mapShape;
 866   TopTools_ListOfShape listShape;
 867
 868   if (aShape.ShapeType() == TopAbs_COMPOUND &&
 869       (TopAbs_ShapeEnum(theShapeType) == TopAbs_SHAPE ||
 870        TopAbs_ShapeEnum(theShapeType) == TopAbs_COMPSOLID ||
 871        TopAbs_ShapeEnum(theShapeType) == TopAbs_COMPOUND)) {
 872     TopoDS_Iterator It (aShape, Standard_True, Standard_True);
 873     for (; It.More(); It.Next()) {
 874       if (mapShape.Add(It.Value())) {
 875         if (TopAbs_ShapeEnum(theShapeType) == TopAbs_SHAPE ||
 876             TopAbs_ShapeEnum(theShapeType) == It.Value().ShapeType()) {
 877           listShape.Append(It.Value());
 878         }
 879       }
 880     }
 881   } else {
 882     TopExp_Explorer exp (aShape, TopAbs_ShapeEnum(theShapeType));
 883     for (; exp.More(); exp.Next())
 884       if (mapShape.Add(exp.Current()))
 885         listShape.Append(exp.Current());
 886   }
 887
 888   if (listShape.IsEmpty()) {
 889     //SetErrorCode("The given shape has no sub-shapes of the requested type");
 890     SetErrorCode(NOT_FOUND_ANY); // NPAL18017
 891     return aSeq;
 892   }
 893
 894   if (isSorted)
 895     SortShapes(listShape);
 896
 897   TopTools_IndexedMapOfShape anIndices;
 898   TopExp::MapShapes(aShape, anIndices);
 899   Handle(TColStd_HArray1OfInteger) anArray;
 900
 901   TopTools_ListIteratorOfListOfShape itSub (listShape);
 902   for (int index = 1; itSub.More(); itSub.Next(), ++index) {
 903     TopoDS_Shape aValue = itSub.Value();
 904     aSeq->Append(anIndices.FindIndex(aValue));
 905   }
 906
 907   Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
 908
 909   //Make a Python command
 910   GEOM::TPythonDump pd (aFunction, /*append=*/true);
 911   pd << "listSubShapeIDs = geompy.SubShapeAll";
 912   pd << (isSorted ? "SortedIDs(" : "IDs(");
 913   pd << theShape << ", " << TopAbs_ShapeEnum(theShapeType) << ")";
 914
 915   SetErrorCode(OK);
 916   return aSeq;
 917 }
 918
 919 //=============================================================================
 920 /*!
 921  *  GetSubShape
 922  */
 923 //=============================================================================
 924 Handle(GEOM_Object) GEOMImpl_IShapesOperations::GetSubShape
 925                                           (Handle(GEOM_Object)    theMainShape,
 926                                            const Standard_Integer theID)
 927 {
 928   SetErrorCode(KO);
 929
 930   if (theMainShape.IsNull()) return NULL;
 931
 932   Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
 933   anArray->SetValue(1, theID);
 934   Handle(GEOM_Object) anObj = GetEngine()->AddSubShape(theMainShape, anArray,true);
 935   if (anObj.IsNull()) {
 936     SetErrorCode("Can not get a sub-shape with the given ID");
 937     return NULL;
 938   }
 939
 940   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
 941
 942   //Make a Python command
 943   GEOM::TPythonDump(aFunction) << anObj << " = geompy.GetSubShape("
 944                                << theMainShape << ", [" << theID << "])";
 945
 946   SetErrorCode(OK);
 947   return anObj;
 948 }
 949
 950 //=============================================================================
 951 /*!
 952  *  GetSubShapeIndex
 953  */
 954 //=============================================================================
 955 Standard_Integer GEOMImpl_IShapesOperations::GetSubShapeIndex (Handle(GEOM_Object) theMainShape,
 956                                                                Handle(GEOM_Object) theSubShape)
 957 {
 958   SetErrorCode(KO);
 959
 960   TopoDS_Shape aMainShape = theMainShape->GetValue();
 961   TopoDS_Shape aSubShape = theSubShape->GetValue();
 962
 963   if (aMainShape.IsNull() || aSubShape.IsNull()) return -1;
 964
 965   TopTools_IndexedMapOfShape anIndices;
 966   TopExp::MapShapes(aMainShape, anIndices);
 967   if (anIndices.Contains(aSubShape)) {
 968     SetErrorCode(OK);
 969     return anIndices.FindIndex(aSubShape);
 970   }
 971
 972   return -1;
 973 }
 974
 975 //=============================================================================
 976 /*!
 977  *  GetTopologyIndex
 978  */
 979 //=============================================================================
 980 Standard_Integer GEOMImpl_IShapesOperations::GetTopologyIndex (Handle(GEOM_Object) theMainShape,
 981                                                                Handle(GEOM_Object) theSubShape)
 982 {
 983   SetErrorCode(OK);
 984
 985   TopoDS_Shape aMainShape = theMainShape->GetValue();
 986   TopoDS_Shape aSubShape = theSubShape->GetValue();
 987
 988   if (aMainShape.IsNull() || aSubShape.IsNull()) {
 989     SetErrorCode("Null argument shape given");
 990     return -1;
 991   }
 992
 993   int index = 1;
 994   if (aSubShape.ShapeType() == TopAbs_COMPOUND) {
 995     TopoDS_Iterator it;
 996     TopTools_ListOfShape CL;
 997     CL.Append(aMainShape);
 998     TopTools_ListIteratorOfListOfShape itC;
 999     for (itC.Initialize(CL); itC.More(); itC.Next()) {
1000       for (it.Initialize(itC.Value()); it.More(); it.Next()) {
1001         if (it.Value().ShapeType() == TopAbs_COMPOUND) {
1002           if (it.Value().IsSame(aSubShape))
1003             return index;
1004           else
1005             index++;
1006           CL.Append(it.Value());
1007         }
1008       }
1009     }
1010   } else {
1011     TopExp_Explorer anExp (aMainShape, aSubShape.ShapeType());
1012     TopTools_MapOfShape M;
1013     for (; anExp.More(); anExp.Next()) {
1014       if (M.Add(anExp.Current())) {
1015         if (anExp.Current().IsSame(aSubShape))
1016           return index;
1017         index++;
1018       }
1019     }
1020   }
1021
1022   SetErrorCode("The sub-shape does not belong to the main shape");
1023   return -1;
1024 }
1025
1026 //=============================================================================
1027 /*!
1028  *  GetShapeTypeString
1029  */
1030 //=============================================================================
1031 TCollection_AsciiString GEOMImpl_IShapesOperations::GetShapeTypeString (Handle(GEOM_Object) theShape)
1032 {
1033   SetErrorCode(KO);
1034
1035   TCollection_AsciiString aTypeName ("Null Shape");
1036
1037   TopoDS_Shape aShape = theShape->GetValue();
1038   if (aShape.IsNull())
1039     return aTypeName;
1040
1041   switch (aShape.ShapeType() )
1042   {
1043   case TopAbs_COMPOUND:
1044     aTypeName = "Compound";
1045     break;
1046   case  TopAbs_COMPSOLID:
1047     aTypeName = "Compound Solid";
1048     break;
1049   case TopAbs_SOLID:
1050     aTypeName = "Solid";
1051     break;
1052   case TopAbs_SHELL:
1053     aTypeName = "Shell";
1054     break;
1055   case TopAbs_FACE:
1056     {
1057       BRepAdaptor_Surface surf (TopoDS::Face(aShape));
1058       if (surf.GetType() == GeomAbs_Plane)
1059         aTypeName = "Plane";
1060       else if (surf.GetType() == GeomAbs_Cylinder)
1061         aTypeName = "Cylindrical Face";
1062       else if (surf.GetType() == GeomAbs_Sphere)
1063         aTypeName = "Spherical Face";
1064       else if (surf.GetType() == GeomAbs_Torus)
1065         aTypeName = "Toroidal Face";
1066       else if (surf.GetType() == GeomAbs_Cone)
1067         aTypeName = "Conical Face";
1068       else
1069         aTypeName = "GEOM::FACE";
1070     }
1071     break;
1072   case TopAbs_WIRE:
1073     aTypeName = "Wire";
1074     break;
1075   case TopAbs_EDGE:
1076     {
1077       BRepAdaptor_Curve curv (TopoDS::Edge(aShape));
1078       if (curv.GetType() == GeomAbs_Line) {
1079         if ((Abs(curv.FirstParameter()) >= 1E6) ||
1080             (Abs(curv.LastParameter()) >= 1E6))
1081           aTypeName = "Line";
1082         else
1083           aTypeName = "Edge" ;
1084       } else if (curv.GetType() == GeomAbs_Circle) {
1085         if (curv.IsClosed())
1086           aTypeName = "Circle";
1087         else
1088           aTypeName = "Arc";
1089       } else {
1090         aTypeName = "Edge";
1091       }
1092     }
1093     break;
1094   case TopAbs_VERTEX:
1095     aTypeName = "Vertex";
1096     break;
1097   case TopAbs_SHAPE:
1098     aTypeName = "Shape";
1099     break;
1100   default:
1101     aTypeName = "Shape of unknown type";
1102   }
1103
1104   return aTypeName;
1105 }
1106
1107
1108 //=============================================================================
1109 /*!
1110  *  NumberOfFaces
1111  */
1112 //=============================================================================
1113 Standard_Integer GEOMImpl_IShapesOperations::NumberOfFaces (Handle(GEOM_Object) theShape)
1114 {
1115   SetErrorCode(KO);
1116
1117   Standard_Integer nb = 0;
1118
1119   if (theShape.IsNull()) return -1;
1120   TopoDS_Shape aShape = theShape->GetValue();
1121   if (aShape.IsNull()) return -1;
1122
1123   TopTools_MapOfShape mapShape;
1124
1125   TopExp_Explorer exp (aShape, TopAbs_FACE);
1126   for (; exp.More(); exp.Next())
1127     if (mapShape.Add(exp.Current()))
1128       nb++;
1129
1130   SetErrorCode(OK);
1131   return nb;
1132 }
1133
1134 //=============================================================================
1135 /*!
1136  *  NumberOfEdges
1137  */
1138 //=============================================================================
1139 Standard_Integer GEOMImpl_IShapesOperations::NumberOfEdges (Handle(GEOM_Object) theShape)
1140 {
1141   SetErrorCode(KO);
1142
1143   Standard_Integer nb = 0;
1144
1145   if (theShape.IsNull()) return -1;
1146   TopoDS_Shape aShape = theShape->GetValue();
1147   if (aShape.IsNull()) return -1;
1148
1149   TopTools_MapOfShape mapShape;
1150
1151   TopExp_Explorer exp (aShape, TopAbs_EDGE);
1152   for (; exp.More(); exp.Next())
1153     if (mapShape.Add(exp.Current()))
1154       nb++;
1155
1156   SetErrorCode(OK);
1157   return nb;
1158 }
1159
1160 //=============================================================================
1161 /*!
1162  *  ReverseShape
1163  */
1164 //=============================================================================
1165 Handle(GEOM_Object) GEOMImpl_IShapesOperations::ReverseShape(Handle(GEOM_Object) theShape)
1166 {
1167   SetErrorCode(KO);
1168
1169   if (theShape.IsNull()) return NULL;
1170
1171   //Add a new reversed object
1172   Handle(GEOM_Object) aReversed = GetEngine()->AddObject(GetDocID(), theShape->GetType());
1173
1174   //Add a new Revese function
1175   Handle(GEOM_Function) aFunction;
1176   aFunction = aReversed->AddFunction(GEOMImpl_ShapeDriver::GetID(), REVERSE_ORIENTATION);
1177   if (aFunction.IsNull()) return NULL;
1178
1179   //Check if the function is set correctly
1180   if (aFunction->GetDriverGUID() != GEOMImpl_ShapeDriver::GetID()) return NULL;
1181
1182   GEOMImpl_IShapes aSI (aFunction);
1183
1184   Handle(GEOM_Function) aRefShape = theShape->GetLastFunction();
1185   if (aRefShape.IsNull()) return NULL;
1186
1187   aSI.SetBase(aRefShape);
1188
1189   //Compute the sub-shape value
1190   try {
1191 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
1192     OCC_CATCH_SIGNALS;
1193 #endif
1194     if (!GetSolver()->ComputeFunction(aFunction)) {
1195       SetErrorCode("Shape driver failed to reverse shape");
1196       return NULL;
1197     }
1198   }
1199   catch (Standard_Failure) {
1200     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1201     SetErrorCode(aFail->GetMessageString());
1202     return NULL;
1203   }
1204
1205   //Make a Python command
1206   GEOM::TPythonDump(aFunction) << aReversed
1207     << " = geompy.ChangeOrientation(" << theShape << ")";
1208
1209   SetErrorCode(OK);
1210   return aReversed;
1211 }
1212
1213 //=============================================================================
1214 /*!
1215  *  GetFreeFacesIDs
1216  */
1217 //=============================================================================
1218 Handle(TColStd_HSequenceOfInteger) GEOMImpl_IShapesOperations::GetFreeFacesIDs
1219                                                  (Handle(GEOM_Object) theShape)
1220 {
1221   SetErrorCode(KO);
1222
1223   if (theShape.IsNull()) return NULL;
1224   TopoDS_Shape aShape = theShape->GetValue();
1225   if (aShape.IsNull()) return NULL;
1226
1227   Handle(TColStd_HSequenceOfInteger) aSeq = new TColStd_HSequenceOfInteger;
1228
1229   TopTools_IndexedDataMapOfShapeListOfShape mapFaceBlocks;
1230   GEOMImpl_Block6Explorer::MapShapesAndAncestors
1231     (aShape, TopAbs_FACE, TopAbs_SOLID, mapFaceBlocks);
1232
1233   Standard_Integer ind = 1, nbFaces = mapFaceBlocks.Extent();
1234
1235   if (nbFaces == 0) {
1236     SetErrorCode("The given shape has no faces");
1237     return aSeq;
1238   }
1239
1240   TopTools_IndexedMapOfShape anIndices;
1241   TopExp::MapShapes(aShape, anIndices);
1242
1243   Standard_Integer id;
1244   for (; ind <= nbFaces; ind++) {
1245     if (mapFaceBlocks.FindFromIndex(ind).Extent() != 2) {
1246       id = anIndices.FindIndex(mapFaceBlocks.FindKey(ind));
1247       aSeq->Append(id);
1248     }
1249   }
1250
1251   //The explode doesn't change object so no new function is required.
1252   Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
1253
1254   //Make a Python command
1255   GEOM::TPythonDump(aFunction, /*append=*/true)
1256     << "listFreeFacesIDs = geompy.GetFreeFacesIDs(" << theShape << ")";
1257
1258   SetErrorCode(OK);
1259   return aSeq;
1260 }
1261
1262 //=======================================================================
1263 //function : GetSharedShapes
1264 //purpose  :
1265 //=======================================================================
1266
1267 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetSharedShapes
1268                                                 (Handle(GEOM_Object)    theShape1,
1269                                                  Handle(GEOM_Object)    theShape2,
1270                                                  const Standard_Integer theShapeType)
1271 {
1272   SetErrorCode(KO);
1273
1274   if (theShape1.IsNull() || theShape2.IsNull()) return NULL;
1275
1276   TopoDS_Shape aShape1 = theShape1->GetValue();
1277   TopoDS_Shape aShape2 = theShape2->GetValue();
1278
1279   if (aShape1.IsNull() || aShape2.IsNull()) return NULL;
1280
1281   TopTools_IndexedMapOfShape anIndices;
1282   TopExp::MapShapes(aShape1, anIndices);
1283   Handle(TColStd_HArray1OfInteger) anArray;
1284
1285   TopTools_IndexedMapOfShape mapShape1;
1286   TopExp::MapShapes(aShape1, TopAbs_ShapeEnum(theShapeType), mapShape1);
1287
1288   Handle(GEOM_Object) anObj;
1289   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
1290   TCollection_AsciiString anAsciiList, anEntry;
1291
1292   TopTools_MapOfShape mapShape2;
1293   TopExp_Explorer exp (aShape2, TopAbs_ShapeEnum(theShapeType));
1294   for (; exp.More(); exp.Next()) {
1295     TopoDS_Shape aSS = exp.Current();
1296     if (mapShape2.Add(aSS) && mapShape1.Contains(aSS)) {
1297       anArray = new TColStd_HArray1OfInteger(1,1);
1298       anArray->SetValue(1, anIndices.FindIndex(aSS));
1299       anObj = GetEngine()->AddSubShape(theShape1, anArray);
1300       aSeq->Append(anObj);
1301
1302       // for python command
1303       TDF_Tool::Entry(anObj->GetEntry(), anEntry);
1304       anAsciiList += anEntry;
1305       anAsciiList += ",";
1306     }
1307   }
1308
1309   if (aSeq->IsEmpty()) {
1310     SetErrorCode("The given shapes have no shared sub-shapes of the requested type");
1311     return aSeq;
1312   }
1313
1314   //Make a Python command
1315   anAsciiList.Trunc(anAsciiList.Length() - 1);
1316
1317   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
1318
1319   GEOM::TPythonDump(aFunction) << "[" << anAsciiList.ToCString()
1320     << "] = geompy.GetSharedShapes(" << theShape1 << ", "
1321       << theShape2 << ", " << TopAbs_ShapeEnum(theShapeType) << ")";
1322
1323   SetErrorCode(OK);
1324   return aSeq;
1325 }
1326
1327 //=============================================================================
1328 /*!
1329  *
1330  */
1331 //=============================================================================
1332 static GEOM::TPythonDump& operator<< (GEOM::TPythonDump&   theDump,
1333                                       const GEOMAlgo_State theState)
1334 {
1335   switch (theState) {
1336   case GEOMAlgo_ST_IN:
1337     theDump << "geompy.GEOM.ST_IN";
1338     break;
1339   case GEOMAlgo_ST_OUT:
1340     theDump << "geompy.GEOM.ST_OUT";
1341     break;
1342   case GEOMAlgo_ST_ON:
1343     theDump << "geompy.GEOM.ST_ON";
1344     break;
1345   case GEOMAlgo_ST_ONIN:
1346     theDump << "geompy.GEOM.ST_ONIN";
1347     break;
1348   case GEOMAlgo_ST_ONOUT:
1349     theDump << "geompy.GEOM.ST_ONOUT";
1350     break;
1351   default:
1352     theDump << "geompy.GEOM.ST_UNKNOWN";
1353     break;
1354   }
1355   return theDump;
1356 }
1357
1358 //=======================================================================
1359 //function : checkTypeShapesOn
1360 /*!
1361  * \brief Checks if theShapeType parameter of GetShapesOnXXX() is OK
1362  * \param theShapeType - the shape type to check
1363  * \retval bool  - result of the check
1364  */
1365 //=======================================================================
1366
1367 bool GEOMImpl_IShapesOperations::checkTypeShapesOn(const Standard_Integer theShapeType)
1368 {
1369   if (theShapeType != TopAbs_VERTEX &&
1370       theShapeType != TopAbs_EDGE &&
1371       theShapeType != TopAbs_FACE &&
1372       theShapeType != TopAbs_SOLID) {
1373     SetErrorCode("Only solids, vertices, edges or faces can be found by this method");
1374     return false;
1375   }
1376   return true;
1377 }
1378
1379 //=======================================================================
1380 //function : makePlane
1381   /*!
1382    * \brief Creates Geom_Plane
1383     * \param theAx1 - shape object defining plane parameters
1384     * \retval Handle(Geom_Surface) - resulting surface
1385    */
1386 //=======================================================================
1387
1388 Handle(Geom_Surface) GEOMImpl_IShapesOperations::makePlane(const TopoDS_Shape& anAx1)
1389 {
1390   if (anAx1.ShapeType() != TopAbs_EDGE) return NULL;
1391   TopoDS_Edge anEdge = TopoDS::Edge(anAx1);
1392   TopoDS_Vertex V1, V2;
1393   TopExp::Vertices(anEdge, V1, V2, Standard_True);
1394   if (V1.IsNull() || V2.IsNull()) {
1395     SetErrorCode("Bad edge given for the plane normal vector");
1396     return NULL;
1397   }
1398   gp_Pnt aLoc = BRep_Tool::Pnt(V1);
1399   gp_Vec aVec (aLoc, BRep_Tool::Pnt(V2));
1400   if (aVec.Magnitude() < Precision::Confusion()) {
1401     SetErrorCode("Vector with null magnitude given");
1402     return NULL;
1403   }
1404   return new Geom_Plane(aLoc, aVec);
1405 }
1406
1407 //=======================================================================
1408 //function : makeCylinder
1409   /*!
1410    * \brief Creates Geom_CylindricalSurface
1411     * \param theAx1 - edge defining cylinder axis
1412     * \param theRadius - cylinder radius
1413     * \retval Handle(Geom_Surface) - resulting surface
1414    */
1415 //=======================================================================
1416
1417 Handle(Geom_Surface) GEOMImpl_IShapesOperations::makeCylinder(const TopoDS_Shape& anAxis,
1418                                                               const Standard_Real theRadius)
1419 {
1420   //Axis of the cylinder
1421   if (anAxis.ShapeType() != TopAbs_EDGE) {
1422     SetErrorCode("Not an edge given for the axis");
1423     return NULL;
1424   }
1425   TopoDS_Edge anEdge = TopoDS::Edge(anAxis);
1426   TopoDS_Vertex V1, V2;
1427   TopExp::Vertices(anEdge, V1, V2, Standard_True);
1428   if (V1.IsNull() || V2.IsNull()) {
1429     SetErrorCode("Bad edge given for the axis");
1430     return NULL;
1431   }
1432   gp_Pnt aLoc = BRep_Tool::Pnt(V1);
1433   gp_Vec aVec (aLoc, BRep_Tool::Pnt(V2));
1434   if (aVec.Magnitude() < Precision::Confusion()) {
1435     SetErrorCode("Vector with null magnitude given");
1436     return NULL;
1437   }
1438
1439   gp_Ax3 anAx3 (aLoc, aVec);
1440   return new Geom_CylindricalSurface(anAx3, theRadius);
1441 }
1442
1443
1444 //=======================================================================
1445 //function : getShapesOnBoxIDs
1446   /*!
1447    * \brief Find IDs of subshapes complying with given status about surface
1448     * \param theBox - the box to check state of subshapes against
1449     * \param theShape - the shape to explore
1450     * \param theShapeType - type of subshape of theShape
1451     * \param theState - required state
1452     * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found subshapes
1453    */
1454 //=======================================================================
1455
1456 Handle(TColStd_HSequenceOfInteger)
1457   GEOMImpl_IShapesOperations::getShapesOnBoxIDs(const Handle(GEOM_Object)& theBox,
1458                                                 const Handle(GEOM_Object)& theShape,
1459                                                 const Standard_Integer theShapeType,
1460                                                 GEOMAlgo_State theState)
1461 {
1462   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs;
1463
1464   TopoDS_Shape aBox = theBox->GetValue();
1465   TopoDS_Shape aShape = theShape->GetValue();
1466
1467   // Check presence of triangulation, build if need
1468   if (!CheckTriangulation(aShape)) {
1469     SetErrorCode("Cannot build triangulation on the shape");
1470     return aSeqOfIDs;
1471   }
1472
1473   // Call algo
1474   GEOMAlgo_FinderShapeOn2 aFinder;
1475   Standard_Real aTol = 0.0001; // default value
1476
1477   Handle(GEOMAlgo_ClsfBox) aClsfBox = new GEOMAlgo_ClsfBox;
1478   aClsfBox->SetBox(aBox);
1479
1480   aFinder.SetShape(aShape);
1481   aFinder.SetTolerance(aTol);
1482   aFinder.SetClsf(aClsfBox);
1483   aFinder.SetShapeType( (TopAbs_ShapeEnum)theShapeType );
1484   aFinder.SetState(theState);
1485   aFinder.Perform();
1486
1487   // Interprete results
1488   Standard_Integer iErr = aFinder.ErrorStatus();
1489   // the detailed description of error codes is in GEOMAlgo_FinderShapeOn1.cxx
1490   if (iErr) {
1491     MESSAGE(" iErr : " << iErr);
1492     TCollection_AsciiString aMsg (" iErr : ");
1493     aMsg += TCollection_AsciiString(iErr);
1494     SetErrorCode(aMsg);
1495     return aSeqOfIDs;
1496   }
1497   Standard_Integer iWrn = aFinder.WarningStatus();
1498   // the detailed description of warning codes is in GEOMAlgo_FinderShapeOn1.cxx
1499   if (iWrn) {
1500     MESSAGE(" *** iWrn : " << iWrn);
1501   }
1502
1503   const TopTools_ListOfShape& listSS = aFinder.Shapes(); // the result
1504
1505   if (listSS.Extent() < 1) {
1506     //SetErrorCode("Not a single sub-shape of the requested type found on the given surface");
1507     SetErrorCode(NOT_FOUND_ANY); // NPAL18017
1508     return aSeqOfIDs;
1509   }
1510
1511   // Fill sequence of object IDs
1512   aSeqOfIDs = new TColStd_HSequenceOfInteger;
1513
1514   TopTools_IndexedMapOfShape anIndices;
1515   TopExp::MapShapes(aShape, anIndices);
1516
1517   TopTools_ListIteratorOfListOfShape itSub (listSS);
1518   for (int index = 1; itSub.More(); itSub.Next(), ++index) {
1519     int id = anIndices.FindIndex(itSub.Value());
1520     aSeqOfIDs->Append(id);
1521   }
1522
1523   return aSeqOfIDs;
1524 }
1525
1526
1527 //=======================================================================
1528 //function : GetShapesOnBoxIDs
1529 /*!
1530    * \brief Find subshapes complying with given status about surface
1531     * \param theBox - the box to check state of subshapes against
1532     * \param theShape - the shape to explore
1533     * \param theShapeType - type of subshape of theShape
1534     * \param theState - required state
1535     * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found subshapes
1536  */
1537 //=======================================================================
1538
1539 Handle(TColStd_HSequenceOfInteger)
1540     GEOMImpl_IShapesOperations::GetShapesOnBoxIDs(const Handle(GEOM_Object)& theBox,
1541                                                   const Handle(GEOM_Object)& theShape,
1542                                                   const Standard_Integer theShapeType,
1543                                                   GEOMAlgo_State theState)
1544 {
1545   // Find subshapes ids
1546   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
1547     getShapesOnBoxIDs (theBox, theShape, theShapeType, theState);
1548   if ( aSeqOfIDs.IsNull()  || aSeqOfIDs->Length() == 0 )
1549     return NULL;
1550
1551   // The GetShapesOnBox() doesn't change object so no new function is required.
1552   Handle(GEOM_Function) aFunction = GEOM::GetCreatedLast(theShape,theBox)->GetLastFunction();
1553
1554   // Make a Python command
1555   GEOM::TPythonDump(aFunction)
1556     << "listShapesOnBoxIDs = geompy.GetShapesOnBoxIDs("
1557     << theBox << ", "
1558     << theShape << ", "
1559     << TopAbs_ShapeEnum(theShapeType) << ", "
1560     << theState << ")";
1561
1562   SetErrorCode(OK);
1563   return aSeqOfIDs;
1564 }
1565
1566 //=======================================================================
1567 //function : GetShapesOnBox
1568 /*!
1569    * \brief Find subshapes complying with given status about surface
1570     * \param theBox - the box to check state of subshapes against
1571     * \param theShape - the shape to explore
1572     * \param theShapeType - type of subshape of theShape
1573     * \param theState - required state
1574     * \retval Handle(TColStd_HSequenceOfTransient) - found subshapes
1575  */
1576 //=======================================================================
1577
1578 Handle(TColStd_HSequenceOfTransient)
1579     GEOMImpl_IShapesOperations::GetShapesOnBox(const Handle(GEOM_Object)& theBox,
1580                                                const Handle(GEOM_Object)&  theShape,
1581                                                const Standard_Integer theShapeType,
1582                                                GEOMAlgo_State theState)
1583 {
1584   // Find subshapes ids
1585   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
1586     getShapesOnBoxIDs (theBox, theShape, theShapeType, theState);
1587   if ( aSeqOfIDs.IsNull()  || aSeqOfIDs->Length() == 0 )
1588     return NULL;
1589
1590   // Find objects by indices
1591   TCollection_AsciiString anAsciiList;
1592   Handle(TColStd_HSequenceOfTransient) aSeq;
1593   aSeq = getObjectsShapesOn( theShape, aSeqOfIDs, anAsciiList );
1594   if ( aSeq.IsNull() || aSeq->IsEmpty() )
1595     return NULL;
1596
1597   // Make a Python command
1598
1599   Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast( aSeq->Value( 1 ));
1600   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
1601
1602   GEOM::TPythonDump(aFunction)
1603     << "[" << anAsciiList.ToCString() << "] = geompy.GetShapesOnBox("
1604     << theBox << ", "
1605     << theShape << ", "
1606     << TopAbs_ShapeEnum(theShapeType) << ", "
1607     << theState << ")";
1608
1609   SetErrorCode(OK);
1610   return aSeq;
1611 }
1612
1613
1614 //=======================================================================
1615 //function : getShapesOnShapeIDs
1616 /*!
1617  * \brief Find IDs of subshapes complying with given status about surface
1618  * \param theCheckShape - the shape to check state of subshapes against
1619  * \param theShape - the shape to explore
1620  * \param theShapeType - type of subshape of theShape
1621  * \param theState - required state
1622  * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found subshapes
1623  */
1624 //=======================================================================
1625
1626 Handle(TColStd_HSequenceOfInteger)
1627   GEOMImpl_IShapesOperations::getShapesOnShapeIDs
1628                                  (const Handle(GEOM_Object)& theCheckShape,
1629                                   const Handle(GEOM_Object)& theShape,
1630                                   const Standard_Integer theShapeType,
1631                                   GEOMAlgo_State theState)
1632 {
1633   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs;
1634
1635   TopoDS_Shape aCheckShape = theCheckShape->GetValue();
1636   TopoDS_Shape aShape = theShape->GetValue();
1637   TopTools_ListOfShape res;
1638
1639   // Check presence of triangulation, build if need
1640   if (!CheckTriangulation(aShape)) {
1641     SetErrorCode("Cannot build triangulation on the shape");
1642     return aSeqOfIDs;
1643   }
1644
1645   // Call algo
1646   GEOMAlgo_FinderShapeOn2 aFinder;
1647   Standard_Real aTol = 0.0001; // default value
1648
1649   Handle(GEOMAlgo_ClsfSolid) aClsfSolid = new GEOMAlgo_ClsfSolid;
1650   aClsfSolid->SetShape(aCheckShape);
1651
1652   aFinder.SetShape(aShape);
1653   aFinder.SetTolerance(aTol);
1654   aFinder.SetClsf(aClsfSolid);
1655   aFinder.SetShapeType( (TopAbs_ShapeEnum)theShapeType );
1656   aFinder.SetState(theState);
1657   aFinder.Perform();
1658
1659   // Interprete results
1660   Standard_Integer iErr = aFinder.ErrorStatus();
1661   // the detailed description of error codes is in GEOMAlgo_FinderShapeOn1.cxx
1662   if (iErr) {
1663     MESSAGE(" iErr : " << iErr);
1664     TCollection_AsciiString aMsg (" iErr : ");
1665     aMsg += TCollection_AsciiString(iErr);
1666     SetErrorCode(aMsg);
1667     return aSeqOfIDs;
1668   }
1669   Standard_Integer iWrn = aFinder.WarningStatus();
1670   // the detailed description of warning codes is in GEOMAlgo_FinderShapeOn1.cxx
1671   if (iWrn) {
1672     MESSAGE(" *** iWrn : " << iWrn);
1673   }
1674
1675   const TopTools_ListOfShape& listSS = aFinder.Shapes(); // the result
1676
1677   if (listSS.Extent() < 1) {
1678     //SetErrorCode("Not a single sub-shape of the requested type found on the given surface");
1679     SetErrorCode(NOT_FOUND_ANY); // NPAL18017
1680   }
1681
1682   // Fill sequence of object IDs
1683   aSeqOfIDs = new TColStd_HSequenceOfInteger;
1684
1685   TopTools_IndexedMapOfShape anIndices;
1686   TopExp::MapShapes(aShape, anIndices);
1687
1688   TopTools_ListIteratorOfListOfShape itSub (listSS);
1689   for (int index = 1; itSub.More(); itSub.Next(), ++index) {
1690     int id = anIndices.FindIndex(itSub.Value());
1691     aSeqOfIDs->Append(id);
1692   }
1693
1694   return aSeqOfIDs;
1695 }
1696
1697
1698 //=======================================================================
1699 //function : GetShapesOnShapeIDs
1700 /*!
1701  * \brief Find subshapes complying with given status about surface
1702  * \param theCheckShape - the shape to check state of subshapes against
1703  * \param theShape - the shape to explore
1704  * \param theShapeType - type of subshape of theShape
1705  * \param theState - required state
1706  * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found subshapes
1707  */
1708 //=======================================================================
1709
1710 Handle(TColStd_HSequenceOfInteger)
1711     GEOMImpl_IShapesOperations::GetShapesOnShapeIDs
1712                             (const Handle(GEOM_Object)& theCheckShape,
1713                              const Handle(GEOM_Object)& theShape,
1714                              const Standard_Integer theShapeType,
1715                              GEOMAlgo_State theState)
1716 {
1717   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
1718     getShapesOnShapeIDs (theCheckShape, theShape, theShapeType, theState);
1719
1720   if ( aSeqOfIDs.IsNull()  || aSeqOfIDs->Length() == 0 )
1721     return NULL;
1722
1723   // The GetShapesOnShape() doesn't change object so no new function is required.
1724   Handle(GEOM_Function) aFunction =
1725     GEOM::GetCreatedLast(theShape,theCheckShape)->GetLastFunction();
1726
1727   // Make a Python command
1728   GEOM::TPythonDump(aFunction)
1729     << "listShapesOnBoxIDs = geompy.GetShapesOnShapeIDs("
1730     << theCheckShape << ", "
1731     << theShape << ", "
1732     << TopAbs_ShapeEnum(theShapeType) << ", "
1733     << theState << ")";
1734
1735   SetErrorCode(OK);
1736   return aSeqOfIDs;
1737 }
1738
1739
1740 //=======================================================================
1741 //function : GetShapesOnShape
1742 /*!
1743  * \brief Find subshapes complying with given status about surface
1744  * \param theCheckShape - the shape to check state of subshapes against
1745  * \param theShape - the shape to explore
1746  * \param theShapeType - type of subshape of theShape
1747  * \param theState - required state
1748  * \retval Handle(TColStd_HSequenceOfTransient) - found subshapes
1749  */
1750 //=======================================================================
1751
1752 Handle(TColStd_HSequenceOfTransient)
1753   GEOMImpl_IShapesOperations::GetShapesOnShape
1754                              (const Handle(GEOM_Object)& theCheckShape,
1755                               const Handle(GEOM_Object)&  theShape,
1756                               const Standard_Integer theShapeType,
1757                               GEOMAlgo_State theState)
1758 {
1759   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
1760     getShapesOnShapeIDs (theCheckShape, theShape, theShapeType, theState);
1761   if ( aSeqOfIDs.IsNull()  || aSeqOfIDs->Length() == 0 )
1762     return NULL;
1763
1764   // Find objects by indices
1765   TCollection_AsciiString anAsciiList;
1766   Handle(TColStd_HSequenceOfTransient) aSeq;
1767   aSeq = getObjectsShapesOn( theShape, aSeqOfIDs, anAsciiList );
1768
1769   if ( aSeq.IsNull() || aSeq->IsEmpty() )
1770     return NULL;
1771
1772   // Make a Python command
1773
1774   Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast( aSeq->Value( 1 ));
1775   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
1776
1777   GEOM::TPythonDump(aFunction)
1778     << "[" << anAsciiList.ToCString() << "] = geompy.GetShapesOnShape("
1779     << theCheckShape << ", "
1780     << theShape << ", "
1781     << TopAbs_ShapeEnum(theShapeType) << ", "
1782     << theState << ")";
1783
1784
1785   SetErrorCode(OK);
1786   return aSeq;
1787 }
1788
1789
1790 //=======================================================================
1791 //function : GetShapesOnShapeAsCompound
1792 //=======================================================================
1793
1794 Handle(GEOM_Object) GEOMImpl_IShapesOperations::GetShapesOnShapeAsCompound
1795                              (const Handle(GEOM_Object)& theCheckShape,
1796                               const Handle(GEOM_Object)&  theShape,
1797                               const Standard_Integer theShapeType,
1798                               GEOMAlgo_State theState)
1799 {
1800   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
1801     getShapesOnShapeIDs (theCheckShape, theShape, theShapeType, theState);
1802
1803   if ( aSeqOfIDs.IsNull()  || aSeqOfIDs->Length() == 0 )
1804     return NULL;
1805
1806   // Find objects by indices
1807   TCollection_AsciiString anAsciiList;
1808   Handle(TColStd_HSequenceOfTransient) aSeq;
1809   aSeq = getObjectsShapesOn( theShape, aSeqOfIDs, anAsciiList );
1810
1811   if ( aSeq.IsNull() || aSeq->IsEmpty() )
1812     return NULL;
1813
1814   TopoDS_Compound aCompound;
1815   BRep_Builder B;
1816   B.MakeCompound(aCompound);
1817   int i = 1;
1818   for(; i<=aSeq->Length(); i++) {
1819     Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast(aSeq->Value(i));
1820     TopoDS_Shape aShape_i = anObj->GetValue();
1821     B.Add(aCompound,aShape_i);
1822   }
1823
1824   //Add a new result object
1825   Handle(GEOM_Object) aRes = GetEngine()->AddObject(GetDocID(), GEOM_SHAPES_ON_SHAPE);
1826   Handle(GEOM_Function) aFunction =
1827     aRes->AddFunction(GEOMImpl_ShapeDriver::GetID(), SHAPES_ON_SHAPE);
1828   aFunction->SetValue(aCompound);
1829
1830   GEOM::TPythonDump(aFunction)
1831     << aRes << " = geompy.GetShapesOnShapeAsCompound("
1832     << theCheckShape << ", "
1833     << theShape << ", "
1834     << TopAbs_ShapeEnum(theShapeType) << ", "
1835     << theState << ")";
1836
1837   SetErrorCode(OK);
1838
1839   return aRes;
1840 }
1841
1842
1843 //=======================================================================
1844 //function : getShapesOnSurfaceIDs
1845   /*!
1846    * \brief Find IDs of subshapes complying with given status about surface
1847     * \param theSurface - the surface to check state of subshapes against
1848     * \param theShape - the shape to explore
1849     * \param theShapeType - type of subshape of theShape
1850     * \param theState - required state
1851     * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found subshapes
1852    */
1853 //=======================================================================
1854
1855 Handle(TColStd_HSequenceOfInteger)
1856   GEOMImpl_IShapesOperations::getShapesOnSurfaceIDs(const Handle(Geom_Surface)& theSurface,
1857                                                     const TopoDS_Shape&         theShape,
1858                                                     TopAbs_ShapeEnum            theShapeType,
1859                                                     GEOMAlgo_State              theState)
1860 {
1861   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs;
1862
1863   // Check presence of triangulation, build if need
1864   if (!CheckTriangulation(theShape)) {
1865     SetErrorCode("Cannot build triangulation on the shape");
1866     return aSeqOfIDs;
1867   }
1868
1869   // Call algo
1870   GEOMAlgo_FinderShapeOn1 aFinder;
1871   Standard_Real aTol = 0.0001; // default value
1872
1873   aFinder.SetShape(theShape);
1874   aFinder.SetTolerance(aTol);
1875   aFinder.SetSurface(theSurface);
1876   aFinder.SetShapeType(theShapeType);
1877   aFinder.SetState(theState);
1878
1879   // Sets the minimal number of inner points for the faces that do not have own
1880   // inner points at all (for e.g. rectangular planar faces have just 2 triangles).
1881   // Default value=3
1882   aFinder.SetNbPntsMin(3);
1883   // Sets the maximal number of inner points for edges or faces.
1884   // It is usefull for the cases when this number is very big (e.g =2000) to improve
1885   // the performance. If this value =0, all inner points will be taken into account.
1886   // Default value=0
1887   aFinder.SetNbPntsMax(100);
1888
1889   aFinder.Perform();
1890
1891   // Interprete results
1892   Standard_Integer iErr = aFinder.ErrorStatus();
1893   // the detailed description of error codes is in GEOMAlgo_FinderShapeOn1.cxx
1894   if (iErr) {
1895     MESSAGE(" iErr : " << iErr);
1896     TCollection_AsciiString aMsg (" iErr : ");
1897     aMsg += TCollection_AsciiString(iErr);
1898     SetErrorCode(aMsg);
1899     return aSeqOfIDs;
1900   }
1901   Standard_Integer iWrn = aFinder.WarningStatus();
1902   // the detailed description of warning codes is in GEOMAlgo_FinderShapeOn1.cxx
1903   if (iWrn) {
1904     MESSAGE(" *** iWrn : " << iWrn);
1905   }
1906
1907   const TopTools_ListOfShape& listSS = aFinder.Shapes(); // the result
1908
1909   if (listSS.Extent() < 1) {
1910     //SetErrorCode("Not a single sub-shape of the requested type found on the given surface");
1911     SetErrorCode(NOT_FOUND_ANY); // NPAL18017
1912     return aSeqOfIDs;
1913   }
1914
1915   // Fill sequence of object IDs
1916   aSeqOfIDs = new TColStd_HSequenceOfInteger;
1917
1918   TopTools_IndexedMapOfShape anIndices;
1919   TopExp::MapShapes(theShape, anIndices);
1920
1921   TopTools_ListIteratorOfListOfShape itSub (listSS);
1922   for (int index = 1; itSub.More(); itSub.Next(), ++index) {
1923     int id = anIndices.FindIndex(itSub.Value());
1924     aSeqOfIDs->Append(id);
1925   }
1926
1927   return aSeqOfIDs;
1928 }
1929
1930 //=======================================================================
1931 //function : getObjectsShapesOn
1932 /*!
1933  * \brief Find shape objects and their entries by their ids
1934  * \param theShapeIDs - incoming shape ids
1935  * \param theShapeEntries - outgoing entries like "entry1, entry2, ..."
1936  * \retval Handle(TColStd_HSequenceOfTransient) - found shape objects
1937  */
1938 //=======================================================================
1939
1940 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::
1941  getObjectsShapesOn(const Handle(GEOM_Object)&                theShape,
1942                     const Handle(TColStd_HSequenceOfInteger)& theShapeIDs,
1943                     TCollection_AsciiString &                 theShapeEntries)
1944 {
1945   Handle(TColStd_HSequenceOfTransient) aSeq;
1946
1947   if ( !theShapeIDs.IsNull() && theShapeIDs->Length() > 0 )
1948   {
1949     aSeq = new TColStd_HSequenceOfTransient;
1950     Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
1951     TCollection_AsciiString anEntry;
1952     for ( int i = 1; i <= theShapeIDs->Length(); ++i )
1953     {
1954       anArray->SetValue(1, theShapeIDs->Value( i ));
1955       Handle(GEOM_Object) anObj = GetEngine()->AddSubShape(theShape, anArray);
1956       aSeq->Append( anObj );
1957
1958       TDF_Tool::Entry(anObj->GetEntry(), anEntry);
1959       if ( i != 1 ) theShapeEntries += ",";
1960       theShapeEntries += anEntry;
1961     }
1962   }
1963   return aSeq;
1964 }
1965
1966 //=======================================================================
1967 //function : getShapesOnSurface
1968 /*!
1969    * \brief Find subshapes complying with given status about surface
1970     * \param theSurface - the surface to check state of subshapes against
1971     * \param theShape - the shape to explore
1972     * \param theShapeType - type of subshape of theShape
1973     * \param theState - required state
1974     * \param theShapeEntries - outgoing entries like "entry1, entry2, ..."
1975     * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found subshapes
1976  */
1977 //=======================================================================
1978
1979 Handle(TColStd_HSequenceOfTransient)
1980     GEOMImpl_IShapesOperations::getShapesOnSurface(const Handle(Geom_Surface)& theSurface,
1981                                                    const Handle(GEOM_Object)&  theShape,
1982                                                    TopAbs_ShapeEnum            theShapeType,
1983                                                    GEOMAlgo_State              theState,
1984                                                    TCollection_AsciiString &   theShapeEntries)
1985 {
1986   // Find subshapes ids
1987   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
1988     getShapesOnSurfaceIDs (theSurface, theShape->GetValue(), theShapeType, theState);
1989   if ( aSeqOfIDs.IsNull()  || aSeqOfIDs->Length() == 0 )
1990     return NULL;
1991
1992   return getObjectsShapesOn( theShape, aSeqOfIDs, theShapeEntries );
1993 }
1994
1995 //=============================================================================
1996 /*!
1997  *  GetShapesOnPlane
1998  */
1999 //=============================================================================
2000 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetShapesOnPlane
2001                                         (const Handle(GEOM_Object)& theShape,
2002                                          const Standard_Integer     theShapeType,
2003                                          const Handle(GEOM_Object)& theAx1,
2004                                          const GEOMAlgo_State       theState)
2005 {
2006   SetErrorCode(KO);
2007
2008   if (theShape.IsNull() || theAx1.IsNull()) return NULL;
2009
2010   TopoDS_Shape aShape = theShape->GetValue();
2011   TopoDS_Shape anAx1  = theAx1->GetValue();
2012
2013   if (aShape.IsNull() || anAx1.IsNull()) return NULL;
2014
2015   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
2016   if ( !checkTypeShapesOn( theShapeType ))
2017     return NULL;
2018
2019   // Create plane
2020   Handle(Geom_Surface) aPlane = makePlane( anAx1 );
2021   if ( aPlane.IsNull() )
2022     return NULL;
2023
2024   // Find objects
2025   TCollection_AsciiString anAsciiList;
2026   Handle(TColStd_HSequenceOfTransient) aSeq;
2027   aSeq = getShapesOnSurface( aPlane, theShape, aShapeType, theState, anAsciiList );
2028   if ( aSeq.IsNull() || aSeq->Length() == 0 )
2029     return NULL;
2030
2031   // Make a Python command
2032
2033   Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast( aSeq->Value( 1 ));
2034   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
2035
2036   GEOM::TPythonDump(aFunction) << "[" << anAsciiList.ToCString()
2037     << "] = geompy.GetShapesOnPlane(" << theShape << ", "
2038       << aShapeType << ", " << theAx1 << ", " << theState << ")";
2039
2040   SetErrorCode(OK);
2041   return aSeq;
2042 }
2043
2044 //=============================================================================
2045 /*!
2046  *  GetShapesOnPlaneWithLocation
2047  */
2048 //=============================================================================
2049 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetShapesOnPlaneWithLocation
2050                                         (const Handle(GEOM_Object)& theShape,
2051                                          const Standard_Integer     theShapeType,
2052                                          const Handle(GEOM_Object)& theAx1,
2053                                          const Handle(GEOM_Object)& thePnt,
2054                                          const GEOMAlgo_State       theState)
2055 {
2056   SetErrorCode(KO);
2057
2058   if (theShape.IsNull() || theAx1.IsNull() || thePnt.IsNull()) return NULL;
2059
2060   TopoDS_Shape aShape = theShape->GetValue();
2061   TopoDS_Shape anAx1  = theAx1->GetValue();
2062   TopoDS_Shape anPnt = thePnt->GetValue();
2063
2064   if (aShape.IsNull() || anAx1.IsNull() || anPnt.IsNull()) return NULL;
2065
2066   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
2067   if ( !checkTypeShapesOn( theShapeType ))
2068     return NULL;
2069
2070   // Create plane
2071   if ( anAx1.ShapeType() != TopAbs_EDGE || anPnt.ShapeType() != TopAbs_VERTEX ) return NULL;
2072   TopoDS_Vertex V1, V2, V3;
2073   TopoDS_Edge anEdge = TopoDS::Edge(anAx1);
2074   TopExp::Vertices(anEdge, V1, V2, Standard_True);
2075
2076   if (V1.IsNull() || V2.IsNull()) {
2077     SetErrorCode("Bad edge given for the plane normal vector");
2078     return NULL;
2079   }
2080   V3 = TopoDS::Vertex(anPnt);
2081
2082   if(V3.IsNull()) {
2083     SetErrorCode("Bad vertex given for the plane location");
2084       return NULL;
2085   }
2086   gp_Pnt aLoc = BRep_Tool::Pnt(V3);
2087   gp_Vec aVec(BRep_Tool::Pnt(V1),BRep_Tool::Pnt(V2));
2088
2089   if (aVec.Magnitude() < Precision::Confusion()) {
2090      SetErrorCode("Vector with null magnitude given");
2091     return NULL;
2092   }
2093   Handle(Geom_Surface) aPlane = new Geom_Plane(aLoc, aVec);
2094
2095   if ( aPlane.IsNull() )
2096     return NULL;
2097
2098   // Find objects
2099   TCollection_AsciiString anAsciiList;
2100   Handle(TColStd_HSequenceOfTransient) aSeq;
2101   aSeq = getShapesOnSurface( aPlane, theShape, aShapeType, theState, anAsciiList );
2102   if ( aSeq.IsNull() || aSeq->Length() == 0 )
2103     return NULL;
2104
2105   // Make a Python command
2106
2107   Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast( aSeq->Value( 1 ));
2108   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
2109
2110   GEOM::TPythonDump(aFunction) << "[" << anAsciiList.ToCString()
2111     << "] = geompy.GetShapesOnPlaneWithLocation(" << theShape << ", "
2112     << aShapeType << ", " << theAx1 << ", "<< thePnt <<", " << theState << ")";
2113
2114   SetErrorCode(OK);
2115   return aSeq;
2116 }
2117
2118 //=============================================================================
2119 /*!
2120  *  GetShapesOnCylinder
2121  */
2122 //=============================================================================
2123 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetShapesOnCylinder
2124                                           (const Handle(GEOM_Object)& theShape,
2125                                            const Standard_Integer     theShapeType,
2126                                            const Handle(GEOM_Object)& theAxis,
2127                                            const Standard_Real        theRadius,
2128                                            const GEOMAlgo_State       theState)
2129 {
2130   SetErrorCode(KO);
2131
2132   if (theShape.IsNull() || theAxis.IsNull()) return NULL;
2133
2134   TopoDS_Shape aShape = theShape->GetValue();
2135   TopoDS_Shape anAxis = theAxis->GetValue();
2136
2137   if (aShape.IsNull() || anAxis.IsNull()) return NULL;
2138
2139   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
2140   if ( !checkTypeShapesOn( aShapeType ))
2141     return NULL;
2142
2143   // Create a cylinder surface
2144   Handle(Geom_Surface) aCylinder = makeCylinder( anAxis, theRadius );
2145   if ( aCylinder.IsNull() )
2146     return NULL;
2147
2148   // Find objects
2149   TCollection_AsciiString anAsciiList;
2150   Handle(TColStd_HSequenceOfTransient) aSeq;
2151   aSeq = getShapesOnSurface( aCylinder, theShape, aShapeType, theState, anAsciiList );
2152   if ( aSeq.IsNull() || aSeq->Length() == 0 )
2153     return NULL;
2154
2155   // Make a Python command
2156
2157   Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast( aSeq->Value( 1 ));
2158   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
2159
2160   GEOM::TPythonDump(aFunction) << "[" << anAsciiList.ToCString()
2161     << "] = geompy.GetShapesOnCylinder(" << theShape << ", " << aShapeType
2162       << ", " << theAxis << ", " << theRadius << ", " << theState << ")";
2163
2164   SetErrorCode(OK);
2165   return aSeq;
2166 }
2167
2168 //=============================================================================
2169 /*!
2170  *  GetShapesOnSphere
2171  */
2172 //=============================================================================
2173 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetShapesOnSphere
2174                                           (const Handle(GEOM_Object)& theShape,
2175                                            const Standard_Integer     theShapeType,
2176                                            const Handle(GEOM_Object)& theCenter,
2177                                            const Standard_Real        theRadius,
2178                                            const GEOMAlgo_State       theState)
2179 {
2180   SetErrorCode(KO);
2181
2182   if (theShape.IsNull() || theCenter.IsNull()) return NULL;
2183
2184   TopoDS_Shape aShape  = theShape->GetValue();
2185   TopoDS_Shape aCenter = theCenter->GetValue();
2186
2187   if (aShape.IsNull() || aCenter.IsNull()) return NULL;
2188
2189   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
2190   if ( !checkTypeShapesOn( aShapeType ))
2191     return NULL;
2192
2193   // Center of the sphere
2194   if (aCenter.ShapeType() != TopAbs_VERTEX) return NULL;
2195   gp_Pnt aLoc = BRep_Tool::Pnt(TopoDS::Vertex(aCenter));
2196
2197   gp_Ax3 anAx3 (aLoc, gp::DZ());
2198   Handle(Geom_SphericalSurface) aSphere =
2199     new Geom_SphericalSurface(anAx3, theRadius);
2200
2201   // Find objects
2202   TCollection_AsciiString anAsciiList;
2203   Handle(TColStd_HSequenceOfTransient) aSeq;
2204   aSeq = getShapesOnSurface( aSphere, theShape, aShapeType, theState, anAsciiList );
2205   if ( aSeq.IsNull() || aSeq->Length() == 0 )
2206     return NULL;
2207
2208   // Make a Python command
2209
2210   Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast( aSeq->Value( 1 ));
2211   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
2212
2213   GEOM::TPythonDump(aFunction) << "[" << anAsciiList.ToCString()
2214     << "] = geompy.GetShapesOnSphere(" << theShape << ", " << aShapeType
2215       << ", " << theCenter << ", " << theRadius << ", " << theState << ")";
2216
2217   SetErrorCode(OK);
2218   return aSeq;
2219 }
2220
2221 //=============================================================================
2222 /*!
2223  *  GetShapesOnPlaneIDs
2224  */
2225 //=============================================================================
2226 Handle(TColStd_HSequenceOfInteger) GEOMImpl_IShapesOperations::GetShapesOnPlaneIDs
2227                                         (const Handle(GEOM_Object)& theShape,
2228                                          const Standard_Integer     theShapeType,
2229                                          const Handle(GEOM_Object)& theAx1,
2230                                          const GEOMAlgo_State       theState)
2231 {
2232   SetErrorCode(KO);
2233
2234   if (theShape.IsNull() || theAx1.IsNull()) return NULL;
2235
2236   TopoDS_Shape aShape = theShape->GetValue();
2237   TopoDS_Shape anAx1  = theAx1->GetValue();
2238
2239   if (aShape.IsNull() || anAx1.IsNull()) return NULL;
2240
2241   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
2242   if ( !checkTypeShapesOn( aShapeType ))
2243     return NULL;
2244
2245   // Create plane
2246   Handle(Geom_Surface) aPlane = makePlane( anAx1 );
2247   if ( aPlane.IsNull() )
2248     return NULL;
2249
2250   // Find object IDs
2251   Handle(TColStd_HSequenceOfInteger) aSeq;
2252   aSeq = getShapesOnSurfaceIDs( aPlane, aShape, aShapeType, theState );
2253
2254   // The GetShapesOnPlaneIDs() doesn't change object so no new function is required.
2255   Handle(GEOM_Function) aFunction = GEOM::GetCreatedLast(theShape,theAx1)->GetLastFunction();
2256
2257   // Make a Python command
2258   GEOM::TPythonDump(aFunction, /*append=*/true)
2259     << "listShapesOnPlane = geompy.GetShapesOnPlaneIDs"
2260     << "(" << theShape << "," << aShapeType << "," << theAx1 << "," << theState << ")";
2261
2262   SetErrorCode(OK);
2263   return aSeq;
2264 }
2265
2266 //=============================================================================
2267 /*!
2268  *  GetShapesOnPlaneWithLocationIDs
2269  */
2270 //=============================================================================
2271 Handle(TColStd_HSequenceOfInteger) GEOMImpl_IShapesOperations::GetShapesOnPlaneWithLocationIDs
2272                                         (const Handle(GEOM_Object)& theShape,
2273                                          const Standard_Integer     theShapeType,
2274                                          const Handle(GEOM_Object)& theAx1,
2275                                          const Handle(GEOM_Object)& thePnt,
2276                                          const GEOMAlgo_State       theState)
2277 {
2278   SetErrorCode(KO);
2279
2280   if (theShape.IsNull() || theAx1.IsNull() || thePnt.IsNull()) return NULL;
2281
2282   TopoDS_Shape aShape = theShape->GetValue();
2283   TopoDS_Shape anAx1  = theAx1->GetValue();
2284   TopoDS_Shape anPnt  = thePnt->GetValue();
2285
2286   if (aShape.IsNull() || anAx1.IsNull() || anPnt.IsNull()) return NULL;
2287
2288   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
2289   if ( !checkTypeShapesOn( aShapeType ))
2290     return NULL;
2291
2292   // Create plane
2293   if (anAx1.ShapeType() != TopAbs_EDGE || anPnt.ShapeType() != TopAbs_VERTEX) return NULL;
2294   TopoDS_Edge anEdge = TopoDS::Edge(anAx1);
2295   TopoDS_Vertex V1, V2, V3;
2296   TopExp::Vertices(anEdge, V1, V2, Standard_True);
2297   if (V1.IsNull() || V2.IsNull()) {
2298     SetErrorCode("Bad edge given for the plane normal vector");
2299     return NULL;
2300   }
2301   V3 = TopoDS::Vertex(anPnt);
2302   if(V3.IsNull()) {
2303     SetErrorCode("Bad vertex given for the plane location");
2304       return NULL;
2305   }
2306   gp_Pnt aLoc = BRep_Tool::Pnt(V3);
2307   gp_Vec aVec(BRep_Tool::Pnt(V1),BRep_Tool::Pnt(V2));
2308   if (aVec.Magnitude() < Precision::Confusion()) {
2309     SetErrorCode("Vector with null magnitude given");
2310     return NULL;
2311   }
2312
2313   Handle(Geom_Surface) aPlane = new Geom_Plane(aLoc, aVec);
2314   if ( aPlane.IsNull() )
2315     return NULL;
2316
2317   // Find object IDs
2318   Handle(TColStd_HSequenceOfInteger) aSeq;
2319   aSeq = getShapesOnSurfaceIDs( aPlane, aShape, aShapeType, theState );
2320
2321   // The GetShapesOnPlaneIDs() doesn't change object so no new function is required.
2322   Handle(GEOM_Function) aFunction = GEOM::GetCreatedLast(theShape,theAx1)->GetLastFunction();
2323
2324   // Make a Python command
2325   GEOM::TPythonDump(aFunction, /*append=*/true)
2326     << "listShapesOnPlane = geompy.GetShapesOnPlaneWithLocationIDs"
2327     << "(" << theShape << ", " << aShapeType << ", " << theAx1 << ", "<< thePnt << ", "  << theState << ")";
2328
2329   SetErrorCode(OK);
2330   return aSeq;
2331 }
2332
2333 //=============================================================================
2334 /*!
2335  *  GetShapesOnCylinderIDs
2336  */
2337 //=============================================================================
2338 Handle(TColStd_HSequenceOfInteger) GEOMImpl_IShapesOperations::GetShapesOnCylinderIDs
2339                                           (const Handle(GEOM_Object)& theShape,
2340                                            const Standard_Integer     theShapeType,
2341                                            const Handle(GEOM_Object)& theAxis,
2342                                            const Standard_Real        theRadius,
2343                                            const GEOMAlgo_State       theState)
2344 {
2345   SetErrorCode(KO);
2346
2347   if (theShape.IsNull() || theAxis.IsNull()) return NULL;
2348
2349   TopoDS_Shape aShape = theShape->GetValue();
2350   TopoDS_Shape anAxis = theAxis->GetValue();
2351
2352   if (aShape.IsNull() || anAxis.IsNull()) return NULL;
2353
2354   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
2355   if ( !checkTypeShapesOn( aShapeType ))
2356     return NULL;
2357
2358   // Create a cylinder surface
2359   Handle(Geom_Surface) aCylinder = makeCylinder( anAxis, theRadius );
2360   if ( aCylinder.IsNull() )
2361     return NULL;
2362
2363   // Find object IDs
2364   Handle(TColStd_HSequenceOfInteger) aSeq;
2365   aSeq = getShapesOnSurfaceIDs( aCylinder, aShape, aShapeType, theState );
2366
2367   // The GetShapesOnCylinder() doesn't change object so no new function is required.
2368   Handle(GEOM_Function) aFunction = GEOM::GetCreatedLast(theShape,theAxis)->GetLastFunction();
2369
2370   // Make a Python command
2371   GEOM::TPythonDump(aFunction, /*append=*/true)
2372     << "listShapesOnCylinder = geompy.GetShapesOnCylinderIDs"
2373     << "(" << theShape << ", " << aShapeType << ", " << theAxis << ", "
2374     << theRadius << ", " << theState << ")";
2375
2376   SetErrorCode(OK);
2377   return aSeq;
2378 }
2379
2380 //=============================================================================
2381 /*!
2382  *  GetShapesOnSphereIDs
2383  */
2384 //=============================================================================
2385 Handle(TColStd_HSequenceOfInteger) GEOMImpl_IShapesOperations::GetShapesOnSphereIDs
2386                                           (const Handle(GEOM_Object)& theShape,
2387                                            const Standard_Integer     theShapeType,
2388                                            const Handle(GEOM_Object)& theCenter,
2389                                            const Standard_Real        theRadius,
2390                                            const GEOMAlgo_State       theState)
2391 {
2392   SetErrorCode(KO);
2393
2394   if (theShape.IsNull() || theCenter.IsNull()) return NULL;
2395
2396   TopoDS_Shape aShape  = theShape->GetValue();
2397   TopoDS_Shape aCenter = theCenter->GetValue();
2398
2399   if (aShape.IsNull() || aCenter.IsNull()) return NULL;
2400
2401   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
2402   if ( !checkTypeShapesOn( aShapeType ))
2403     return NULL;
2404
2405   // Center of the sphere
2406   if (aCenter.ShapeType() != TopAbs_VERTEX) return NULL;
2407   gp_Pnt aLoc = BRep_Tool::Pnt(TopoDS::Vertex(aCenter));
2408
2409   gp_Ax3 anAx3 (aLoc, gp::DZ());
2410   Handle(Geom_SphericalSurface) aSphere =
2411     new Geom_SphericalSurface(anAx3, theRadius);
2412
2413   // Find object IDs
2414   Handle(TColStd_HSequenceOfInteger) aSeq;
2415   aSeq = getShapesOnSurfaceIDs( aSphere, aShape, aShapeType, theState );
2416
2417   // The GetShapesOnSphere() doesn't change object so no new function is required.
2418   Handle(GEOM_Function) aFunction = GEOM::GetCreatedLast(theShape,theCenter)->GetLastFunction();
2419
2420   // Make a Python command
2421   GEOM::TPythonDump(aFunction, /*append=*/true)
2422     << "listShapesOnCylinder = geompy.GetShapesOnCylinderIDs"
2423     << "(" << theShape << ", " << aShapeType << ", " << theCenter << ", "
2424     << theRadius << ", " << theState << ")";
2425
2426   SetErrorCode(OK);
2427   return aSeq;
2428 }
2429
2430 //=======================================================================
2431 //function : getShapesOnQuadrangleIDs
2432   /*!
2433    * \brief Find IDs of subshapes complying with given status about quadrangle
2434     * \param theShape - the shape to explore
2435     * \param theShapeType - type of subshape of theShape
2436     * \param theTopLeftPoint - top left quadrangle corner
2437     * \param theTopRigthPoint - top right quadrangle corner
2438     * \param theBottomLeftPoint - bottom left quadrangle corner
2439     * \param theBottomRigthPoint - bottom right quadrangle corner
2440     * \param theState - required state
2441     * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found subshapes
2442    */
2443 //=======================================================================
2444
2445 Handle(TColStd_HSequenceOfInteger)
2446   GEOMImpl_IShapesOperations::getShapesOnQuadrangleIDs (const Handle(GEOM_Object)& theShape,
2447                                                         const Standard_Integer     theShapeType,
2448                                                         const Handle(GEOM_Object)& theTopLeftPoint,
2449                                                         const Handle(GEOM_Object)& theTopRigthPoint,
2450                                                         const Handle(GEOM_Object)& theBottomLeftPoint,
2451                                                         const Handle(GEOM_Object)& theBottomRigthPoint,
2452                                                         const GEOMAlgo_State       theState)
2453 {
2454   SetErrorCode(KO);
2455
2456   if ( theShape.IsNull() ||
2457        theTopLeftPoint.IsNull() ||
2458        theTopRigthPoint.IsNull() ||
2459        theBottomLeftPoint.IsNull() ||
2460        theBottomRigthPoint.IsNull() )
2461     return NULL;
2462
2463   TopoDS_Shape aShape = theShape->GetValue();
2464   TopoDS_Shape aTL = theTopLeftPoint->GetValue();
2465   TopoDS_Shape aTR = theTopRigthPoint->GetValue();
2466   TopoDS_Shape aBL = theBottomLeftPoint->GetValue();
2467   TopoDS_Shape aBR = theBottomRigthPoint->GetValue();
2468
2469   if (aShape.IsNull() ||
2470       aTL.IsNull() ||
2471       aTR.IsNull() ||
2472       aBL.IsNull() ||
2473       aBR.IsNull() ||
2474       aTL.ShapeType() != TopAbs_VERTEX ||
2475       aTR.ShapeType() != TopAbs_VERTEX ||
2476       aBL.ShapeType() != TopAbs_VERTEX ||
2477       aBR.ShapeType() != TopAbs_VERTEX )
2478     return NULL;
2479
2480   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
2481   if ( !checkTypeShapesOn( aShapeType ))
2482     return NULL;
2483
2484   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs;
2485
2486   // Check presence of triangulation, build if need
2487   if (!CheckTriangulation(aShape)) {
2488     SetErrorCode("Cannot build triangulation on the shape");
2489     return aSeqOfIDs;
2490   }
2491
2492   // Call algo
2493   gp_Pnt aPntTL = BRep_Tool::Pnt(TopoDS::Vertex(aTL));
2494   gp_Pnt aPntTR = BRep_Tool::Pnt(TopoDS::Vertex(aTR));
2495   gp_Pnt aPntBL = BRep_Tool::Pnt(TopoDS::Vertex(aBL));
2496   gp_Pnt aPntBR = BRep_Tool::Pnt(TopoDS::Vertex(aBR));
2497
2498   GEOMAlgo_FinderShapeOnQuad aFinder( aPntTL, aPntTR, aPntBL, aPntBR );
2499   Standard_Real aTol = 0.0001; // default value
2500
2501   aFinder.SetShape(aShape);
2502   aFinder.SetTolerance(aTol);
2503   //aFinder.SetSurface(theSurface);
2504   aFinder.SetShapeType(aShapeType);
2505   aFinder.SetState(theState);
2506
2507   // Sets the minimal number of inner points for the faces that do not have own
2508   // inner points at all (for e.g. rectangular planar faces have just 2 triangles).
2509   // Default value=3
2510   aFinder.SetNbPntsMin(3);
2511   // Sets the maximal number of inner points for edges or faces.
2512   // It is usefull for the cases when this number is very big (e.g =2000) to improve
2513   // the performance. If this value =0, all inner points will be taken into account.
2514   // Default value=0
2515   aFinder.SetNbPntsMax(100);
2516
2517   aFinder.Perform();
2518
2519   // Interprete results
2520   Standard_Integer iErr = aFinder.ErrorStatus();
2521   // the detailed description of error codes is in GEOMAlgo_FinderShapeOn1.cxx
2522   if (iErr) {
2523     MESSAGE(" iErr : " << iErr);
2524     TCollection_AsciiString aMsg (" iErr : ");
2525     aMsg += TCollection_AsciiString(iErr);
2526     SetErrorCode(aMsg);
2527     return aSeqOfIDs;
2528   }
2529   Standard_Integer iWrn = aFinder.WarningStatus();
2530   // the detailed description of warning codes is in GEOMAlgo_FinderShapeOn1.cxx
2531   if (iWrn) {
2532     MESSAGE(" *** iWrn : " << iWrn);
2533   }
2534
2535   const TopTools_ListOfShape& listSS = aFinder.Shapes(); // the result
2536
2537   if (listSS.Extent() < 1) {
2538     //SetErrorCode("Not a single sub-shape of the requested type found on the given surface");
2539     SetErrorCode(NOT_FOUND_ANY); // NPAL18017
2540     return aSeqOfIDs;
2541   }
2542
2543   // Fill sequence of object IDs
2544   aSeqOfIDs = new TColStd_HSequenceOfInteger;
2545
2546   TopTools_IndexedMapOfShape anIndices;
2547   TopExp::MapShapes(aShape, anIndices);
2548
2549   TopTools_ListIteratorOfListOfShape itSub (listSS);
2550   for (int index = 1; itSub.More(); itSub.Next(), ++index) {
2551     int id = anIndices.FindIndex(itSub.Value());
2552     aSeqOfIDs->Append(id);
2553   }
2554   return aSeqOfIDs;
2555 }
2556
2557 //=======================================================================
2558 //function : GetShapesOnQuadrangle
2559   /*!
2560    * \brief Find subshapes complying with given status about quadrangle
2561     * \param theShape - the shape to explore
2562     * \param theShapeType - type of subshape of theShape
2563     * \param theTopLeftPoint - top left quadrangle corner
2564     * \param theTopRigthPoint - top right quadrangle corner
2565     * \param theBottomLeftPoint - bottom left quadrangle corner
2566     * \param theBottomRigthPoint - bottom right quadrangle corner
2567     * \param theState - required state
2568     * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found subshapes
2569    */
2570 //=======================================================================
2571
2572 Handle(TColStd_HSequenceOfTransient)
2573     GEOMImpl_IShapesOperations::GetShapesOnQuadrangle (const Handle(GEOM_Object)& theShape,
2574                                                        const Standard_Integer     theShapeType,
2575                                                        const Handle(GEOM_Object)& theTopLeftPoint,
2576                                                        const Handle(GEOM_Object)& theTopRigthPoint,
2577                                                        const Handle(GEOM_Object)& theBottomLeftPoint,
2578                                                        const Handle(GEOM_Object)& theBottomRigthPoint,
2579                                                        const GEOMAlgo_State       theState)
2580 {
2581   // Find indices
2582   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
2583     getShapesOnQuadrangleIDs( theShape,
2584                               theShapeType,
2585                               theTopLeftPoint,
2586                               theTopRigthPoint,
2587                               theBottomLeftPoint,
2588                               theBottomRigthPoint,
2589                               theState);
2590   if ( aSeqOfIDs.IsNull() || aSeqOfIDs->IsEmpty() )
2591     return NULL;
2592
2593   // Find objects by indices
2594   TCollection_AsciiString anAsciiList;
2595   Handle(TColStd_HSequenceOfTransient) aSeq;
2596   aSeq = getObjectsShapesOn( theShape, aSeqOfIDs, anAsciiList );
2597   if ( aSeq.IsNull() || aSeq->IsEmpty() )
2598     return NULL;
2599
2600   // Make a Python command
2601
2602   Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast( aSeq->Value( 1 ));
2603   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
2604
2605   GEOM::TPythonDump(aFunction)
2606     << "[" << anAsciiList.ToCString() << "] = geompy.GetShapesOnQuadrangle("
2607     << theShape << ", "
2608     << TopAbs_ShapeEnum(theShapeType) << ", "
2609     << theTopLeftPoint << ", "
2610     << theTopRigthPoint << ", "
2611     << theBottomLeftPoint << ", "
2612     << theBottomRigthPoint << ", "
2613     << theState << ")";
2614
2615   SetErrorCode(OK);
2616   return aSeq;
2617 }
2618
2619 //=======================================================================
2620 //function : GetShapesOnQuadrangleIDs
2621   /*!
2622    * \brief Find IDs of subshapes complying with given status about quadrangle
2623     * \param theShape - the shape to explore
2624     * \param theShapeType - type of subshape of theShape
2625     * \param theTopLeftPoint - top left quadrangle corner
2626     * \param theTopRigthPoint - top right quadrangle corner
2627     * \param theBottomLeftPoint - bottom left quadrangle corner
2628     * \param theBottomRigthPoint - bottom right quadrangle corner
2629     * \param theState - required state
2630     * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found subshapes
2631    */
2632 //=======================================================================
2633
2634 Handle(TColStd_HSequenceOfInteger)
2635   GEOMImpl_IShapesOperations::GetShapesOnQuadrangleIDs (const Handle(GEOM_Object)& theShape,
2636                                                         const Standard_Integer     theShapeType,
2637                                                         const Handle(GEOM_Object)& theTopLeftPoint,
2638                                                         const Handle(GEOM_Object)& theTopRigthPoint,
2639                                                         const Handle(GEOM_Object)& theBottomLeftPoint,
2640                                                         const Handle(GEOM_Object)& theBottomRigthPoint,
2641                                                         const GEOMAlgo_State       theState)
2642 {
2643   // Find indices
2644   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
2645     getShapesOnQuadrangleIDs( theShape,
2646                               theShapeType,
2647                               theTopLeftPoint,
2648                               theTopRigthPoint,
2649                               theBottomLeftPoint,
2650                               theBottomRigthPoint,
2651                               theState);
2652   if ( aSeqOfIDs.IsNull() || aSeqOfIDs->IsEmpty() )
2653     return NULL;
2654
2655   // Make a Python command
2656
2657   // The GetShapesOnCylinder() doesn't change object so no new function is required.
2658   Handle(GEOM_Object) lastObj = GEOM::GetCreatedLast(theShape,theTopLeftPoint);
2659   lastObj = GEOM::GetCreatedLast(lastObj,theTopRigthPoint);
2660   lastObj = GEOM::GetCreatedLast(lastObj,theBottomRigthPoint);
2661   lastObj = GEOM::GetCreatedLast(lastObj,theBottomLeftPoint);
2662   Handle(GEOM_Function) aFunction = lastObj->GetLastFunction();
2663
2664   GEOM::TPythonDump(aFunction, /*append=*/true)
2665     << "listShapesOnQuadrangle = geompy.GetShapesOnQuadrangleIDs("
2666     << theShape << ", "
2667     << TopAbs_ShapeEnum(theShapeType) << ", "
2668     << theTopLeftPoint << ", "
2669     << theTopRigthPoint << ", "
2670     << theBottomLeftPoint << ", "
2671     << theBottomRigthPoint << ", "
2672     << theState << ")";
2673
2674   SetErrorCode(OK);
2675   return aSeqOfIDs;
2676 }
2677
2678
2679 //=============================================================================
2680 /*!
2681  *  GetInPlaceOfShape
2682  */
2683 //=============================================================================
2684 static bool GetInPlaceOfShape (const Handle(GEOM_Function)& theWhereFunction,
2685                                const TopTools_IndexedMapOfShape& theWhereIndices,
2686                                const TopoDS_Shape& theWhat,
2687                                TColStd_ListOfInteger& theModifiedList)
2688 {
2689   if (theWhereFunction.IsNull() || theWhat.IsNull()) return false;
2690
2691   if (theWhereIndices.Contains(theWhat)) {
2692     // entity was not changed by the operation
2693     Standard_Integer aWhatIndex = theWhereIndices.FindIndex(theWhat);
2694     theModifiedList.Append(aWhatIndex);
2695     return true;
2696   }
2697
2698   // try to find in history
2699   TDF_Label aHistoryLabel = theWhereFunction->GetHistoryEntry(Standard_False);
2700
2701   // search in history for all argument shapes
2702   Standard_Boolean isFound = Standard_False;
2703   Standard_Boolean isGood = Standard_False;
2704
2705   TDF_LabelSequence aLabelSeq;
2706   theWhereFunction->GetDependency(aLabelSeq);
2707   Standard_Integer nbArg = aLabelSeq.Length();
2708
2709   for (Standard_Integer iarg = 1; iarg <= nbArg && !isFound; iarg++) {
2710
2711     TDF_Label anArgumentRefLabel = aLabelSeq.Value(iarg);
2712
2713     Handle(GEOM_Object) anArgumentObject = GEOM_Object::GetReferencedObject(anArgumentRefLabel);
2714     TopoDS_Shape anArgumentShape = anArgumentObject->GetValue();
2715
2716     TopTools_IndexedMapOfShape anArgumentIndices;
2717     TopExp::MapShapes(anArgumentShape, anArgumentIndices);
2718
2719     if (anArgumentIndices.Contains(theWhat)) {
2720       isFound = Standard_True;
2721       Standard_Integer aWhatIndex = anArgumentIndices.FindIndex(theWhat);
2722
2723       // Find corresponding label in history
2724       TDF_Label anArgumentHistoryLabel =
2725         theWhereFunction->GetArgumentHistoryEntry(anArgumentRefLabel, Standard_False);
2726       if (anArgumentHistoryLabel.IsNull()) {
2727         // Lost History of operation argument. Possibly, all its entities was removed.
2728         isGood = Standard_True;
2729       }
2730       else {
2731         TDF_Label aWhatHistoryLabel = anArgumentHistoryLabel.FindChild(aWhatIndex, Standard_False);
2732
2733         if (aWhatHistoryLabel.IsNull()) {
2734           // Removed entity ? Compound ? Compsolid ? Shell ? Wire
2735           isGood = Standard_False;
2736         } else {
2737           Handle(TDataStd_IntegerArray) anIntegerArray;
2738           if (!aWhatHistoryLabel.FindAttribute(TDataStd_IntegerArray::GetID(), anIntegerArray)) {
2739             //Error: Empty modifications history for the sought shape.
2740             isGood = Standard_False;
2741           }
2742           else {
2743             isGood = Standard_True;
2744             Standard_Integer imod, aModifLen = anIntegerArray->Array()->Length();
2745             for (imod = 1; imod <= aModifLen; imod++) {
2746               theModifiedList.Append(anIntegerArray->Array()->Value(imod));
2747             }
2748           }
2749         }
2750       }
2751     }
2752   }
2753
2754   isFound = isGood;
2755
2756   if (!isFound) {
2757     // try compound/compsolid/shell/wire element by element
2758     bool isFoundAny = false;
2759     TopTools_MapOfShape mapShape;
2760
2761     if (theWhat.ShapeType() == TopAbs_COMPOUND ||
2762         theWhat.ShapeType() == TopAbs_COMPSOLID) {
2763       // recursive processing of compound/compsolid
2764       TopoDS_Iterator anIt (theWhat, Standard_True, Standard_True);
2765       for (; anIt.More(); anIt.Next()) {
2766         if (mapShape.Add(anIt.Value())) {
2767           TopoDS_Shape curWhat = anIt.Value();
2768           isFoundAny = GetInPlaceOfShape(theWhereFunction, theWhereIndices, curWhat, theModifiedList);
2769           if (isFoundAny) isFound = Standard_True;
2770         }
2771       }
2772     }
2773     else if (theWhat.ShapeType() == TopAbs_SHELL) {
2774       // try to replace a shell by its faces images
2775       TopExp_Explorer anExp (theWhat, TopAbs_FACE);
2776       for (; anExp.More(); anExp.Next()) {
2777         if (mapShape.Add(anExp.Current())) {
2778           TopoDS_Shape curWhat = anExp.Current();
2779           isFoundAny = GetInPlaceOfShape(theWhereFunction, theWhereIndices, curWhat, theModifiedList);
2780           if (isFoundAny) isFound = Standard_True;
2781         }
2782       }
2783     }
2784     else if (theWhat.ShapeType() == TopAbs_WIRE) {
2785       // try to replace a wire by its edges images
2786       TopExp_Explorer anExp (theWhat, TopAbs_EDGE);
2787       for (; anExp.More(); anExp.Next()) {
2788         if (mapShape.Add(anExp.Current())) {
2789           TopoDS_Shape curWhat = anExp.Current();
2790           isFoundAny = GetInPlaceOfShape(theWhereFunction, theWhereIndices, curWhat, theModifiedList);
2791           if (isFoundAny) isFound = Standard_True;
2792         }
2793       }
2794     }
2795     else {
2796       // Removed entity
2797     }
2798   }
2799
2800   return isFound;
2801 }
2802
2803 //=============================================================================
2804 /*!
2805  *  GetShapeProperties
2806  */
2807 //=============================================================================
2808
2809 void GEOMImpl_IShapesOperations::GetShapeProperties( const TopoDS_Shape aShape, Standard_Real tab[],
2810                                                      gp_Pnt & aVertex )
2811 {
2812   GProp_GProps theProps;
2813   gp_Pnt aCenterMass;
2814   TopoDS_Shape aPntShape;
2815   Standard_Real aShapeSize;
2816
2817   if      (aShape.ShapeType() == TopAbs_EDGE) BRepGProp::LinearProperties(aShape,  theProps);
2818   else if (aShape.ShapeType() == TopAbs_FACE) BRepGProp::SurfaceProperties(aShape, theProps);
2819   else                                        BRepGProp::VolumeProperties(aShape,  theProps);
2820
2821   aCenterMass = theProps.CentreOfMass();
2822   aShapeSize  = theProps.Mass();
2823
2824   aPntShape = BRepBuilderAPI_MakeVertex(aCenterMass).Shape();
2825   aVertex   = BRep_Tool::Pnt( TopoDS::Vertex( aPntShape ) );
2826   tab[0] = aVertex.X();
2827   tab[1] = aVertex.Y();
2828   tab[2] = aVertex.Z();
2829   tab[3] = aShapeSize;
2830   return;
2831 }
2832
2833 //=============================================================================
2834 /*!
2835  *  GetInPlace
2836  */
2837 //=============================================================================
2838 Handle(GEOM_Object) GEOMImpl_IShapesOperations::GetInPlace (Handle(GEOM_Object) theShapeWhere,
2839                                                             Handle(GEOM_Object) theShapeWhat)
2840 {
2841   SetErrorCode(KO);
2842
2843   if (theShapeWhere.IsNull() || theShapeWhat.IsNull()) return NULL;
2844
2845   TopoDS_Shape aWhere = theShapeWhere->GetValue();
2846   TopoDS_Shape aWhat  = theShapeWhat->GetValue();
2847   TopoDS_Shape aPntShape;
2848   TopoDS_Vertex aVertex;
2849
2850   if (aWhere.IsNull() || aWhat.IsNull()) {
2851     SetErrorCode("Error: aWhere and aWhat TopoDS_Shape are Null.");
2852     return NULL;
2853   }
2854
2855   Handle(GEOM_Function) aWhereFunction = theShapeWhere->GetLastFunction();
2856   if (aWhereFunction.IsNull()) {
2857     SetErrorCode("Error: aWhereFunction is Null.");
2858     return NULL;
2859   }
2860
2861   TopTools_IndexedMapOfShape aWhereIndices;
2862   TopExp::MapShapes(aWhere, aWhereIndices);
2863
2864   TColStd_ListOfInteger aModifiedList;
2865   Standard_Integer aWhereIndex;
2866   Handle(TColStd_HArray1OfInteger) aModifiedArray;
2867   Handle(GEOM_Object) aResult;
2868
2869   bool isFound = false;
2870   Standard_Integer iType = TopAbs_SOLID;
2871   Standard_Integer compType = TopAbs_SOLID;
2872   Standard_Real    aWhat_Mass = 0., aWhere_Mass = 0.;
2873   Standard_Real    tab_aWhat[4],    tab_aWhere[4];
2874   Standard_Real    dl_l = 1e-3;
2875   Standard_Real    min_l, Tol_1D, Tol_2D, Tol_3D, Tol_Mass;
2876   Standard_Real    aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
2877   Bnd_Box          BoundingBox;
2878   gp_Pnt           aPnt, aPnt_aWhat, tab_Pnt[2];
2879   GProp_GProps     aProps;
2880
2881   // Find the iType of the aWhat shape
2882   if      ( aWhat.ShapeType() == TopAbs_EDGE  || aWhat.ShapeType() == TopAbs_WIRE )      iType = TopAbs_EDGE;
2883   else if ( aWhat.ShapeType() == TopAbs_FACE  || aWhat.ShapeType() == TopAbs_SHELL )     iType = TopAbs_FACE;
2884   else if ( aWhat.ShapeType() == TopAbs_SOLID || aWhat.ShapeType() == TopAbs_COMPSOLID ) iType = TopAbs_SOLID;
2885   else if ( aWhat.ShapeType() == TopAbs_COMPOUND ) {
2886     // Only the iType of the first shape in the compound is taken into account
2887     TopoDS_Iterator It (aWhat, Standard_True, Standard_True);
2888     compType = It.Value().ShapeType();
2889     if      ( compType == TopAbs_EDGE  || compType == TopAbs_WIRE )     iType = TopAbs_EDGE;
2890     else if ( compType == TopAbs_FACE  || compType == TopAbs_SHELL)     iType = TopAbs_FACE;
2891     else if ( compType == TopAbs_SOLID || compType == TopAbs_COMPSOLID) iType = TopAbs_SOLID;
2892   }
2893   else {
2894     SetErrorCode("Error: An attempt to extract a shape of not supported type.");
2895     return NULL;
2896   }
2897
2898   TopExp_Explorer Exp_aWhat( aWhat,   TopAbs_ShapeEnum( iType ) );
2899   TopExp_Explorer Exp_aWhere( aWhere, TopAbs_ShapeEnum( iType ) );
2900   TopExp_Explorer Exp_Edge( aWhere,   TopAbs_EDGE );
2901
2902   // Find the shortest edge in theShapeWhere shape
2903   BRepBndLib::Add(aWhere, BoundingBox);
2904   BoundingBox.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
2905   min_l = fabs(aXmax - aXmin);
2906   if( min_l < fabs(aYmax - aYmin) ) min_l = fabs(aYmax - aYmin);
2907   if( min_l < fabs(aZmax - aZmin) ) min_l = fabs(aZmax - aZmin);
2908   min_l /= dl_l;
2909   for ( Standard_Integer nbEdge = 0; Exp_Edge.More(); Exp_Edge.Next(), nbEdge++ ) {
2910     TopExp_Explorer Exp_Vertex( Exp_Edge.Current(), TopAbs_VERTEX);
2911     for ( Standard_Integer nbVertex = 0; Exp_Vertex.More(); Exp_Vertex.Next(), nbVertex++ ) {
2912       aPnt = BRep_Tool::Pnt( TopoDS::Vertex( Exp_Vertex.Current() ) );
2913       tab_Pnt[nbVertex] = aPnt;
2914     }
2915     if ( ! tab_Pnt[0].IsEqual(tab_Pnt[1], dl_l) ) {
2916       BRepGProp::LinearProperties(Exp_Edge.Current(), aProps);
2917       if ( aProps.Mass() < min_l ) min_l = aProps.Mass();
2918     }
2919   }
2920
2921   // Compute tolerances
2922   Tol_1D = dl_l * min_l;
2923   Tol_2D = dl_l * ( min_l * min_l) * ( 2. + dl_l);
2924   Tol_3D = dl_l * ( min_l * min_l * min_l ) * ( 3. + (3 * dl_l) + (dl_l * dl_l) );
2925
2926   Tol_Mass = Tol_3D;
2927   if      ( iType == TopAbs_EDGE ) Tol_Mass = Tol_1D;
2928   else if ( iType == TopAbs_FACE ) Tol_Mass = Tol_2D;
2929
2930   // Compute the ShapeWhat Mass
2931   for ( ; Exp_aWhat.More(); Exp_aWhat.Next() ) {
2932     if      ( iType == TopAbs_EDGE ) BRepGProp::LinearProperties(Exp_aWhat.Current(),  aProps);
2933     else if ( iType == TopAbs_FACE ) BRepGProp::SurfaceProperties(Exp_aWhat.Current(), aProps);
2934     else                             BRepGProp::VolumeProperties(Exp_aWhat.Current(),  aProps);
2935     aWhat_Mass += aProps.Mass();
2936   }
2937
2938   // Searching for the sub-shapes inside the ShapeWhere shape
2939   for ( Exp_aWhere.ReInit(); Exp_aWhere.More(); Exp_aWhere.Next() ) {
2940     GetShapeProperties( Exp_aWhere.Current(), tab_aWhere, aPnt );
2941     for ( Exp_aWhat.ReInit(); Exp_aWhat.More(); Exp_aWhat.Next() ) {
2942       GetShapeProperties( Exp_aWhat.Current(), tab_aWhat, aPnt_aWhat );
2943       if ( fabs(tab_aWhat[3] - tab_aWhere[3]) <= Tol_Mass && aPnt_aWhat.Distance(aPnt) <= Tol_1D ) isFound = true;
2944       else {
2945         if ( (tab_aWhat[3] - tab_aWhere[3]) > Tol_Mass ) {
2946           aPntShape = BRepBuilderAPI_MakeVertex( aPnt ).Shape();
2947           aVertex   = TopoDS::Vertex( aPntShape );
2948           BRepExtrema_DistShapeShape aWhereDistance ( aVertex, Exp_aWhere.Current() );
2949           BRepExtrema_DistShapeShape aWhatDistance  ( aVertex, Exp_aWhat.Current() );
2950           if ( fabs(aWhereDistance.Value() - aWhatDistance.Value()) <= Tol_1D ) isFound = true;
2951         }
2952       }
2953       if ( isFound ) {
2954         aWhereIndex = aWhereIndices.FindIndex(Exp_aWhere.Current());
2955         aModifiedList.Append(aWhereIndex);
2956         aWhere_Mass += tab_aWhere[3];
2957         isFound = false;
2958         break;
2959       }
2960     }
2961     if ( fabs( aWhat_Mass - aWhere_Mass ) <= Tol_Mass ) break;
2962   }
2963
2964   if (aModifiedList.Extent() == 0) { // Not found any Results
2965     SetErrorCode(NOT_FOUND_ANY);
2966     return NULL;
2967   }
2968
2969   aModifiedArray = new TColStd_HArray1OfInteger (1, aModifiedList.Extent());
2970   TColStd_ListIteratorOfListOfInteger anIterModif (aModifiedList);
2971   for (Standard_Integer imod = 1; anIterModif.More(); anIterModif.Next(), imod++)
2972     aModifiedArray->SetValue(imod, anIterModif.Value());
2973
2974   //Add a new object
2975   aResult = GetEngine()->AddSubShape(theShapeWhere, aModifiedArray);
2976   if (aResult.IsNull()) {
2977     SetErrorCode("Error in algorithm: result found, but cannot be returned.");
2978     return NULL;
2979   }
2980
2981   if (aModifiedArray->Length() > 1) {
2982     //Set a GROUP type
2983     aResult->SetType(GEOM_GROUP);
2984
2985     //Set a sub shape type
2986     TopoDS_Shape aFirstFound = aWhereIndices.FindKey(aModifiedArray->Value(1));
2987     TopAbs_ShapeEnum aShapeType = aFirstFound.ShapeType();
2988
2989     TDF_Label aFreeLabel = aResult->GetFreeLabel();
2990     TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)aShapeType);
2991   }
2992
2993   //Make a Python command
2994   Handle(GEOM_Function) aFunction = aResult->GetFunction(1);
2995
2996   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetInPlace("
2997     << theShapeWhere << ", " << theShapeWhat << ")";
2998
2999   SetErrorCode(OK);
3000   return aResult;
3001 }
3002
3003 //=======================================================================
3004 //function : GetInPlaceByHistory
3005 //purpose  :
3006 //=======================================================================
3007 Handle(GEOM_Object) GEOMImpl_IShapesOperations::GetInPlaceByHistory
3008                                           (Handle(GEOM_Object) theShapeWhere,
3009                                            Handle(GEOM_Object) theShapeWhat)
3010 {
3011   SetErrorCode(KO);
3012
3013   if (theShapeWhere.IsNull() || theShapeWhat.IsNull()) return NULL;
3014
3015   TopoDS_Shape aWhere = theShapeWhere->GetValue();
3016   TopoDS_Shape aWhat  = theShapeWhat->GetValue();
3017
3018   if (aWhere.IsNull() || aWhat.IsNull()) return NULL;
3019
3020   Handle(GEOM_Function) aWhereFunction = theShapeWhere->GetLastFunction();
3021   if (aWhereFunction.IsNull()) return NULL;
3022
3023   //Fill array of indices
3024   TopTools_IndexedMapOfShape aWhereIndices;
3025   TopExp::MapShapes(aWhere, aWhereIndices);
3026
3027   // process shape
3028   TColStd_ListOfInteger aModifiedList;
3029   bool isFound = GetInPlaceOfShape(aWhereFunction, aWhereIndices, aWhat, aModifiedList);
3030
3031   if (!isFound || aModifiedList.Extent() < 1) {
3032     SetErrorCode("Error: No history found for the sought shape or its sub-shapes.");
3033     return NULL;
3034   }
3035
3036   Handle(TColStd_HArray1OfInteger) aModifiedArray =
3037     new TColStd_HArray1OfInteger (1, aModifiedList.Extent());
3038   TColStd_ListIteratorOfListOfInteger anIterModif (aModifiedList);
3039   for (Standard_Integer imod = 1; anIterModif.More(); anIterModif.Next(), imod++) {
3040     aModifiedArray->SetValue(imod, anIterModif.Value());
3041   }
3042
3043   //Add a new object
3044   Handle(GEOM_Object) aResult = GetEngine()->AddSubShape(theShapeWhere, aModifiedArray);
3045   if (aResult.IsNull()) {
3046     SetErrorCode("Error in algorithm: result found, but cannot be returned.");
3047     return NULL;
3048   }
3049
3050   if (aModifiedArray->Length() > 1) {
3051     //Set a GROUP type
3052     aResult->SetType(GEOM_GROUP);
3053
3054     //Set a sub shape type
3055     TopoDS_Shape aFirstFound = aWhereIndices.FindKey(aModifiedArray->Value(1));
3056     TopAbs_ShapeEnum aShapeType = aFirstFound.ShapeType();
3057
3058     TDF_Label aFreeLabel = aResult->GetFreeLabel();
3059     TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)aShapeType);
3060   }
3061
3062   //Make a Python command
3063   Handle(GEOM_Function) aFunction = aResult->GetFunction(1);
3064
3065   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetInPlace("
3066     << theShapeWhere << ", " << theShapeWhat << ")";
3067
3068   SetErrorCode(OK);
3069   return aResult;
3070 }
3071
3072 //=======================================================================
3073 //function : SortShapes
3074 //purpose  :
3075 //=======================================================================
3076 void GEOMImpl_IShapesOperations::SortShapes(TopTools_ListOfShape& SL)
3077 {
3078   Standard_Integer MaxShapes = SL.Extent();
3079   TopTools_Array1OfShape  aShapes (1,MaxShapes);
3080   TColStd_Array1OfInteger OrderInd(1,MaxShapes);
3081   TColStd_Array1OfReal    MidXYZ  (1,MaxShapes); //X,Y,Z;
3082   TColStd_Array1OfReal    Length  (1,MaxShapes); //X,Y,Z;
3083
3084   // Computing of CentreOfMass
3085   Standard_Integer Index;
3086   GProp_GProps GPr;
3087   gp_Pnt GPoint;
3088   TopTools_ListIteratorOfListOfShape it(SL);
3089   for (Index=1;  it.More();  Index++)
3090   {
3091     TopoDS_Shape S = it.Value();
3092     SL.Remove( it ); // == it.Next()
3093     aShapes(Index) = S;
3094     OrderInd.SetValue (Index, Index);
3095     if (S.ShapeType() == TopAbs_VERTEX)
3096     {
3097       GPoint = BRep_Tool::Pnt( TopoDS::Vertex( S ));
3098       Length.SetValue( Index, (Standard_Real) S.Orientation());
3099     }
3100     else
3101     {
3102       BRepGProp::LinearProperties (S, GPr);
3103       GPoint = GPr.CentreOfMass();
3104       Length.SetValue( Index, GPr.Mass() );
3105     }
3106     MidXYZ.SetValue(Index,
3107                     GPoint.X()*999 + GPoint.Y()*99 + GPoint.Z()*0.9);
3108     //cout << Index << " L: " << Length(Index) << "CG: " << MidXYZ(Index) << endl;
3109   }
3110
3111   // Sorting
3112   Standard_Integer aTemp;
3113   Standard_Boolean exchange, Sort = Standard_True;
3114   Standard_Real    tol = Precision::Confusion();
3115   while (Sort)
3116   {
3117     Sort = Standard_False;
3118     for (Index=1; Index < MaxShapes; Index++)
3119     {
3120       Standard_Real dMidXYZ = MidXYZ(OrderInd(Index)) - MidXYZ(OrderInd(Index+1));
3121       Standard_Real dLength = Length(OrderInd(Index)) - Length(OrderInd(Index+1));
3122       if ( dMidXYZ >= tol ) {
3123 //         cout << "MidXYZ: " << MidXYZ(OrderInd(Index))<< " > " <<MidXYZ(OrderInd(Index+1))
3124 //              << " d: " << dMidXYZ << endl;
3125         exchange = Standard_True;
3126       }
3127       else if ( Abs(dMidXYZ) < tol && dLength >= tol ) {
3128 //         cout << "Length: " << Length(OrderInd(Index))<< " > " <<Length(OrderInd(Index+1))
3129 //              << " d: " << dLength << endl;
3130         exchange = Standard_True;
3131       }
3132       else if ( Abs(dMidXYZ) < tol && Abs(dLength) < tol &&
3133                 aShapes(OrderInd(Index)).ShapeType() <= TopAbs_FACE) {
3134         // PAL17233
3135         // equal values possible on shapes such as two halves of a sphere and
3136         // a membrane inside the sphere
3137         Bnd_Box box1,box2;
3138         BRepBndLib::Add( aShapes( OrderInd(Index) ), box1 );
3139         if ( box1.IsVoid() ) continue;
3140         BRepBndLib::Add( aShapes( OrderInd(Index+1) ), box2 );
3141         Standard_Real dSquareExtent = box1.SquareExtent() - box2.SquareExtent();
3142         if ( dSquareExtent >= tol ) {
3143 //           cout << "SquareExtent: " << box1.SquareExtent()<<" > "<<box2.SquareExtent() << endl;
3144           exchange = Standard_True;
3145         }
3146         else if ( Abs(dSquareExtent) < tol ) {
3147           Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, val1, val2;
3148           box1.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
3149           val1 = (aXmin+aXmax)*999 + (aYmin+aYmax)*99 + (aZmin+aZmax)*0.9;
3150           box2.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
3151           val2 = (aXmin+aXmax)*999 + (aYmin+aYmax)*99 + (aZmin+aZmax)*0.9;
3152           exchange = val1 > val2;
3153 //           cout << "box: " << val1<<" > "<<val2 << endl;
3154         }
3155       }
3156       else {
3157         exchange = Standard_False;
3158       }
3159       if (exchange)
3160       {
3161 //         cout << "exchange " << Index << " & " << Index+1 << endl;
3162         aTemp = OrderInd(Index);
3163         OrderInd(Index) = OrderInd(Index+1);
3164         OrderInd(Index+1) = aTemp;
3165         Sort = Standard_True;
3166       }
3167     }
3168   }
3169
3170   for (Index=1; Index <= MaxShapes; Index++)
3171     SL.Append( aShapes( OrderInd(Index) ));
3172 }
3173
3174 //=======================================================================
3175 //function : CompsolidToCompound
3176 //purpose  :
3177 //=======================================================================
3178 TopoDS_Shape GEOMImpl_IShapesOperations::CompsolidToCompound (const TopoDS_Shape& theCompsolid)
3179 {
3180   if (theCompsolid.ShapeType() != TopAbs_COMPSOLID) {
3181     return theCompsolid;
3182   }
3183
3184   TopoDS_Compound aCompound;
3185   BRep_Builder B;
3186   B.MakeCompound(aCompound);
3187
3188   TopTools_MapOfShape mapShape;
3189   TopoDS_Iterator It (theCompsolid, Standard_True, Standard_True);
3190
3191   for (; It.More(); It.Next()) {
3192     TopoDS_Shape aShape_i = It.Value();
3193     if (mapShape.Add(aShape_i)) {
3194       B.Add(aCompound, aShape_i);
3195     }
3196   }
3197
3198   return aCompound;
3199 }
3200
3201 //=======================================================================
3202 //function : CheckTriangulation
3203 //purpose  :
3204 //=======================================================================
3205 bool GEOMImpl_IShapesOperations::CheckTriangulation (const TopoDS_Shape& aShape)
3206 {
3207   bool isTriangulation = true;
3208
3209   TopExp_Explorer exp (aShape, TopAbs_FACE);
3210   if (exp.More())
3211   {
3212     TopLoc_Location aTopLoc;
3213     Handle(Poly_Triangulation) aTRF;
3214     aTRF = BRep_Tool::Triangulation(TopoDS::Face(exp.Current()), aTopLoc);
3215     if (aTRF.IsNull()) {
3216       isTriangulation = false;
3217     }
3218   }
3219   else // no faces, try edges
3220   {
3221     TopExp_Explorer expe (aShape, TopAbs_EDGE);
3222     if (!expe.More()) {
3223       return false;
3224     }
3225     TopLoc_Location aLoc;
3226     Handle(Poly_Polygon3D) aPE = BRep_Tool::Polygon3D(TopoDS::Edge(expe.Current()), aLoc);
3227     if (aPE.IsNull()) {
3228       isTriangulation = false;
3229     }
3230   }
3231
3232   if (!isTriangulation) {
3233     // calculate deflection
3234     Standard_Real aDeviationCoefficient = 0.001;
3235
3236     Bnd_Box B;
3237     BRepBndLib::Add(aShape, B);
3238     Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
3239     B.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
3240
3241     Standard_Real dx = aXmax - aXmin, dy = aYmax - aYmin, dz = aZmax - aZmin;
3242     Standard_Real aDeflection = Max(Max(dx, dy), dz) * aDeviationCoefficient * 4;
3243     Standard_Real aHLRAngle = 0.349066;
3244
3245     BRepMesh_IncrementalMesh Inc (aShape, aDeflection, Standard_False, aHLRAngle);
3246   }
3247
3248   return true;
3249 }
3250
3251 #define MAX_TOLERANCE 1.e-7
3252
3253
3254 //=======================================================================
3255 //function : isSameEdge
3256 //purpose  : Returns True if two edges coincide
3257 //=======================================================================
3258 static bool isSameEdge(const TopoDS_Edge& theEdge1, const TopoDS_Edge& theEdge2)
3259 {
3260   TopoDS_Vertex V11, V12, V21, V22;
3261   TopExp::Vertices(theEdge1, V11, V12);
3262   TopExp::Vertices(theEdge2, V21, V22);
3263   gp_Pnt P11 = BRep_Tool::Pnt(V11);
3264   gp_Pnt P12 = BRep_Tool::Pnt(V12);
3265   gp_Pnt P21 = BRep_Tool::Pnt(V21);
3266   gp_Pnt P22 = BRep_Tool::Pnt(V22);
3267   bool coincide = false;
3268
3269   //Check that ends of edges coincide
3270   if(P11.Distance(P21) <= MAX_TOLERANCE) {
3271     if(P12.Distance(P22) <= MAX_TOLERANCE) coincide =  true;
3272   }
3273   else if(P11.Distance(P22) <= MAX_TOLERANCE) {
3274     if(P12.Distance(P21) <= MAX_TOLERANCE) coincide = true;
3275   }
3276
3277   if(!coincide) return false;
3278
3279   if (BRep_Tool::Degenerated(theEdge1))
3280     if (BRep_Tool::Degenerated(theEdge2)) return true;
3281     else return false;
3282   else
3283     if (BRep_Tool::Degenerated(theEdge2)) return false;
3284
3285   double U11, U12, U21, U22;
3286   Handle(Geom_Curve) C1 = BRep_Tool::Curve(theEdge1, U11, U12);
3287   Handle(Geom_Curve) C2 = BRep_Tool::Curve(theEdge2, U21, U22);
3288   if(C1->DynamicType() == C2->DynamicType()) return true;
3289
3290   //Check that both edges has the same geometry
3291   double range = U12-U11;
3292   double U = U11+ range/3.0;
3293   gp_Pnt P1 = C1->Value(U);     //Compute a point on one third of the edge's length
3294   U = U11+range*2.0/3.0;
3295   gp_Pnt P2 = C1->Value(U);     //Compute a point on two thirds of the edge's length
3296
3297   if(!GeomLib_Tool::Parameter(C2, P1, MAX_TOLERANCE, U) ||  U < U21 || U > U22)
3298     return false;
3299
3300   if(P1.Distance(C2->Value(U)) > MAX_TOLERANCE) return false;
3301
3302   if(!GeomLib_Tool::Parameter(C2, P2, MAX_TOLERANCE, U) || U < U21 || U > U22)
3303     return false;
3304
3305   if(P2.Distance(C2->Value(U)) > MAX_TOLERANCE) return false;
3306
3307   return true;
3308 }
3309
3310 #include <TopoDS_TShape.hxx>
3311 //=======================================================================
3312 //function : isSameFace
3313 //purpose  : Returns True if two faces coincide
3314 //=======================================================================
3315 static bool isSameFace(const TopoDS_Face& theFace1, const TopoDS_Face& theFace2)
3316 {
3317   TopExp_Explorer E(theFace1, TopAbs_EDGE);
3318   TopTools_ListOfShape LS1, LS2;
3319   for(; E.More(); E.Next()) LS1.Append(E.Current());
3320
3321   E.Init(theFace2, TopAbs_EDGE);
3322   for(; E.More(); E.Next()) LS2.Append(E.Current());
3323
3324   //Compare the number of edges in the faces
3325   if(LS1.Extent() != LS2.Extent()) return false;
3326
3327   double aMin = RealFirst(), aMax = RealLast();
3328   double xminB1=aMax, yminB1=aMax, zminB1=aMax, xminB2=aMax, yminB2=aMax, zminB2=aMax;
3329   double xmaxB1=aMin, ymaxB1=aMin, zmaxB1=aMin, xmaxB2=aMin, ymaxB2=aMin, zmaxB2=aMin;
3330
3331   for(E.Init(theFace1, TopAbs_VERTEX); E.More(); E.Next()) {
3332     gp_Pnt P = BRep_Tool::Pnt(TopoDS::Vertex(E.Current()));
3333     if(P.X() < xminB1) xminB1 = P.X();
3334     if(P.Y() < yminB1) yminB1 = P.Y();
3335     if(P.Z() < zminB1) zminB1 = P.Z();
3336     if(P.X() > xmaxB1) xmaxB1 = P.X();
3337     if(P.Y() > ymaxB1) ymaxB1 = P.Y();
3338     if(P.Z() > zmaxB1) zmaxB1 = P.Z();
3339   }
3340
3341   for(E.Init(theFace2, TopAbs_VERTEX); E.More(); E.Next()) {
3342     gp_Pnt P = BRep_Tool::Pnt(TopoDS::Vertex(E.Current()));
3343     if(P.X() < xminB2) xminB2 = P.X();
3344     if(P.Y() < yminB2) yminB2 = P.Y();
3345     if(P.Z() < zminB2) zminB2 = P.Z();
3346     if(P.X() > xmaxB2) xmaxB2 = P.X();
3347     if(P.Y() > ymaxB2) ymaxB2 = P.Y();
3348     if(P.Z() > zmaxB2) zmaxB2 = P.Z();
3349   }
3350
3351
3352   //Compare the bounding boxes of both faces
3353   if(gp_Pnt(xminB1, yminB1, zminB1).Distance(gp_Pnt(xminB2, yminB2, zminB2)) > MAX_TOLERANCE)
3354     return false;
3355
3356   if(gp_Pnt(xmaxB1, ymaxB1, zmaxB1).Distance(gp_Pnt(xmaxB2, ymaxB2, zmaxB2)) > MAX_TOLERANCE)
3357     return false;
3358
3359   Handle(Geom_Surface) S1 = BRep_Tool::Surface(theFace1);
3360   Handle(Geom_Surface) S2 = BRep_Tool::Surface(theFace2);
3361
3362   //Check if there a coincidence of two surfaces at least in two points
3363   double U11, U12, V11, V12, U21, U22, V21, V22;
3364   BRepTools::UVBounds(theFace1, U11, U12, V11, V12);
3365   BRepTools::UVBounds(theFace2, U21, U22, V21, V22);
3366
3367   double rangeU = U12-U11;
3368   double rangeV = V12-V11;
3369   double U = U11 + rangeU/3.0;
3370   double V = V11 + rangeV/3.0;
3371   gp_Pnt P1 = S1->Value(U, V);
3372   U = U11+rangeU*2.0/3.0;
3373   V = V11+rangeV*2.0/3.0;
3374   gp_Pnt P2 = S1->Value(U, V);
3375
3376   if(!GeomLib_Tool::Parameters(S2, P1, MAX_TOLERANCE, U, V) || U < U21 || U > U22 || V < V21 || V > V22)
3377     return false;
3378
3379   if(P1.Distance(S2->Value(U,V)) > MAX_TOLERANCE) return false;
3380
3381   if(!GeomLib_Tool::Parameters(S2, P2, MAX_TOLERANCE, U, V) || U < U21 || U > U22 || V < V21 || V > V22)
3382     return false;
3383
3384   if(P2.Distance(S2->Value(U, V)) > MAX_TOLERANCE) return false;
3385
3386   //Check that each edge of the Face1 has a counterpart in the Face2
3387   TopTools_MapOfOrientedShape aMap;
3388   TopTools_ListIteratorOfListOfShape LSI1(LS1);
3389   for(; LSI1.More(); LSI1.Next()) {
3390     TopoDS_Edge E = TopoDS::Edge(LSI1.Value());
3391     bool isFound = false;
3392     TopTools_ListIteratorOfListOfShape LSI2(LS2);
3393     for(; LSI2.More(); LSI2.Next()) {
3394       TopoDS_Shape aValue = LSI2.Value();
3395       if(aMap.Contains(aValue)) continue; //To avoid checking already found edge several times
3396       if(isSameEdge(E, TopoDS::Edge(aValue))) {
3397         aMap.Add(aValue);
3398         isFound = true;
3399         break;
3400       }
3401     }
3402     if(!isFound) return false;
3403   }
3404
3405   return true;
3406 }
3407
3408 //=======================================================================
3409 //function : isSameSolid
3410 //purpose  : Returns True if two solids coincide
3411 //=======================================================================
3412 bool isSameSolid(const TopoDS_Solid& theSolid1, const TopoDS_Solid& theSolid2)
3413 {
3414   TopExp_Explorer E(theSolid1, TopAbs_FACE);
3415   TopTools_ListOfShape LS1, LS2;
3416   for(; E.More(); E.Next()) LS1.Append(E.Current());
3417   E.Init(theSolid2, TopAbs_FACE);
3418   for(; E.More(); E.Next()) LS2.Append(E.Current());
3419
3420   if(LS1.Extent() != LS2.Extent()) return false;
3421
3422   double aMin = RealFirst(), aMax = RealLast();
3423   double xminB1=aMax, yminB1=aMax, zminB1=aMax, xminB2=aMax, yminB2=aMax, zminB2=aMax;
3424   double xmaxB1=aMin, ymaxB1=aMin, zmaxB1=aMin, xmaxB2=aMin, ymaxB2=aMin, zmaxB2=aMin;
3425
3426   for(E.Init(theSolid1, TopAbs_VERTEX); E.More(); E.Next()) {
3427     gp_Pnt P = BRep_Tool::Pnt(TopoDS::Vertex(E.Current()));
3428     if(P.X() < xminB1) xminB1 = P.X();
3429     if(P.Y() < yminB1) yminB1 = P.Y();
3430     if(P.Z() < zminB1) zminB1 = P.Z();
3431     if(P.X() > xmaxB1) xmaxB1 = P.X();
3432     if(P.Y() > ymaxB1) ymaxB1 = P.Y();
3433     if(P.Z() > zmaxB1) zmaxB1 = P.Z();
3434   }
3435
3436   for(E.Init(theSolid2, TopAbs_VERTEX); E.More(); E.Next()) {
3437     gp_Pnt P = BRep_Tool::Pnt(TopoDS::Vertex(E.Current()));
3438     if(P.X() < xminB2) xminB2 = P.X();
3439     if(P.Y() < yminB2) yminB2 = P.Y();
3440     if(P.Z() < zminB2) zminB2 = P.Z();
3441     if(P.X() > xmaxB2) xmaxB2 = P.X();
3442     if(P.Y() > ymaxB2) ymaxB2 = P.Y();
3443     if(P.Z() > zmaxB2) zmaxB2 = P.Z();
3444   }
3445
3446   //Compare the bounding boxes of both solids
3447   if(gp_Pnt(xminB1, yminB1, zminB1).Distance(gp_Pnt(xminB2, yminB2, zminB2)) > MAX_TOLERANCE)
3448     return false;
3449
3450   if(gp_Pnt(xmaxB1, ymaxB1, zmaxB1).Distance(gp_Pnt(xmaxB2, ymaxB2, zmaxB2)) > MAX_TOLERANCE)
3451     return false;
3452
3453   //Check that each face of the Solid1 has a counterpart in the Solid2
3454   TopTools_MapOfOrientedShape aMap;
3455   TopTools_ListIteratorOfListOfShape LSI1(LS1);
3456   for(; LSI1.More(); LSI1.Next()) {
3457     TopoDS_Face F = TopoDS::Face(LSI1.Value());
3458     bool isFound = false;
3459     TopTools_ListIteratorOfListOfShape LSI2(LS2);
3460     for(; LSI2.More(); LSI2.Next()) {
3461       if(aMap.Contains(LSI2.Value())) continue; //To avoid checking already found faces several times
3462       if(isSameFace(F, TopoDS::Face(LSI2.Value()))) {
3463         aMap.Add(LSI2.Value());
3464         isFound = true;
3465         break;
3466       }
3467     }
3468     if(!isFound) return false;
3469   }
3470
3471   return true;
3472 }
3473
3474 //=======================================================================
3475 //function : GetSame
3476 //purpose  :
3477 //=======================================================================
3478 Handle(GEOM_Object) GEOMImpl_IShapesOperations::GetSame(const Handle(GEOM_Object)& theShapeWhere,
3479                                                         const Handle(GEOM_Object)& theShapeWhat)
3480 {
3481   SetErrorCode(KO);
3482   if (theShapeWhere.IsNull() || theShapeWhat.IsNull()) return NULL;
3483
3484   TopoDS_Shape aWhere = theShapeWhere->GetValue();
3485   TopoDS_Shape aWhat  = theShapeWhat->GetValue();
3486
3487   if (aWhere.IsNull() || aWhat.IsNull()) return NULL;
3488
3489   int anIndex = -1;
3490   bool isFound = false;
3491   TopoDS_Shape aSubShape;
3492   TopTools_MapOfShape aMap;
3493
3494   switch(aWhat.ShapeType()) {
3495     case TopAbs_VERTEX: {
3496       gp_Pnt P = BRep_Tool::Pnt(TopoDS::Vertex(aWhat));
3497       TopExp_Explorer E(aWhere, TopAbs_VERTEX);
3498       for(; E.More(); E.Next()) {
3499         if(!aMap.Add(E.Current())) continue;
3500         gp_Pnt P2 = BRep_Tool::Pnt(TopoDS::Vertex(E.Current()));
3501         if(P.Distance(P2) <= MAX_TOLERANCE) {
3502           isFound = true;
3503           aSubShape = E.Current();
3504           break;
3505         }
3506       }
3507       break;
3508                         }
3509     case TopAbs_FACE: {
3510       TopoDS_Face aFace = TopoDS::Face(aWhat);
3511       TopExp_Explorer E(aWhere, TopAbs_FACE);
3512       for(; E.More(); E.Next()) {
3513         if(!aMap.Add(E.Current())) continue;
3514         if(isSameFace(aFace, TopoDS::Face(E.Current()))) {
3515           aSubShape = E.Current();
3516           isFound = true;
3517           break;
3518         }
3519       }
3520       break;
3521                       }
3522     case TopAbs_EDGE: {
3523       TopoDS_Edge anEdge = TopoDS::Edge(aWhat);
3524       TopExp_Explorer E(aWhere, TopAbs_EDGE);
3525       for(; E.More(); E.Next()) {
3526         if(!aMap.Add(E.Current())) continue;
3527         if(isSameEdge(anEdge, TopoDS::Edge(E.Current()))) {
3528           aSubShape = E.Current();
3529           isFound = true;
3530           break;
3531         }
3532       }
3533       break;
3534                       }
3535     case TopAbs_SOLID: {
3536       TopoDS_Solid aSolid = TopoDS::Solid(aWhat);
3537       TopExp_Explorer E(aWhere, TopAbs_SOLID);
3538       for(; E.More(); E.Next()) {
3539         if(!aMap.Add(E.Current())) continue;
3540         if(isSameSolid(aSolid, TopoDS::Solid(E.Current()))) {
3541           aSubShape = E.Current();
3542           isFound = true;
3543           break;
3544         }
3545       }
3546       break;
3547                        }
3548     default:
3549       return NULL;
3550   }
3551
3552   if(isFound) {
3553     TopTools_IndexedMapOfShape anIndices;
3554     TopExp::MapShapes(aWhere, anIndices);
3555     if (anIndices.Contains(aSubShape))
3556       anIndex = anIndices.FindIndex(aSubShape);
3557   }
3558
3559   if(anIndex < 0) return NULL;
3560
3561   Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
3562
3563   anArray->SetValue(1, anIndex);
3564
3565   Handle(GEOM_Object) aResult = GetEngine()->AddSubShape(theShapeWhere, anArray);
3566   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
3567
3568   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetSame("
3569     << theShapeWhere << ", " << theShapeWhat << ")";
3570
3571   SetErrorCode(OK);
3572
3573   return aResult;
3574 }