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