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