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