src/GEOMImpl/GEOMImpl_IShapesOperations.cxx

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