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