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