src/GEOMImpl/GEOMImpl_IShapesOperations.cxx

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