src/GEOMImpl/GEOMImpl_IShapesOperations.cxx

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