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