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