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