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