src/GEOMImpl/GEOMImpl_IShapesOperations.cxx

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