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