Salome HOME
b81b261ce0e0a7c62376ca50427225407460d986
[modules/geom.git] / src / GEOMImpl / GEOMImpl_IShapesOperations.cxx
1 // Copyright (C) 2007-2012  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
23 //  File      : GEOMImpl_IShapesOperations.cxx
24 //  Created   :
25 //  Author    : modified by Lioka RAZAFINDRAZAKA (CEA) 22/06/2007
26 //  Project   : SALOME
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 #include "GEOMImpl_IHealingOperations.hxx"
45
46 #include <GEOMImpl_Gen.hxx>
47
48 #include "GEOM_Function.hxx"
49 #include "GEOM_ISubShape.hxx"
50 #include "GEOM_PythonDump.hxx"
51
52 #include "GEOMAlgo_ClsfBox.hxx"
53 #include "GEOMAlgo_ClsfSolid.hxx"
54 #include "GEOMAlgo_CoupleOfShapes.hxx"
55 #include "GEOMAlgo_FinderShapeOn1.hxx"
56 #include "GEOMAlgo_FinderShapeOnQuad.hxx"
57 #include "GEOMAlgo_FinderShapeOn2.hxx"
58 #include "GEOMAlgo_GetInPlace.hxx"
59 #include "GEOMAlgo_GlueDetector.hxx"
60 #include "GEOMAlgo_ListIteratorOfListOfCoupleOfShapes.hxx"
61 #include "GEOMAlgo_ListOfCoupleOfShapes.hxx"
62
63 #include <Basics_OCCTVersion.hxx>
64
65 #include "utilities.h"
66 #include "OpUtil.hxx"
67 #include "Utils_ExceptHandlers.hxx"
68
69 #include <TFunction_DriverTable.hxx>
70 #include <TFunction_Driver.hxx>
71 #include <TFunction_Logbook.hxx>
72 #include <TDataStd_Integer.hxx>
73 #include <TDataStd_IntegerArray.hxx>
74 #include <TDataStd_ListIteratorOfListOfExtendedString.hxx>
75 #include <TDF_Tool.hxx>
76
77 #include <BRepExtrema_ExtCF.hxx>
78 #include <BRepExtrema_DistShapeShape.hxx>
79
80 #include <BRep_Tool.hxx>
81 #include <BRep_Builder.hxx>
82 #include <BRepTools.hxx>
83 #include <BRepGProp.hxx>
84 #include <BRepAdaptor_Curve.hxx>
85 #include <BRepAdaptor_Surface.hxx>
86 #include <BRepBndLib.hxx>
87 #include <BRepMesh_IncrementalMesh.hxx>
88
89 #include <TopAbs.hxx>
90 #include <TopExp.hxx>
91 #include <TopExp_Explorer.hxx>
92 #include <TopLoc_Location.hxx>
93 #include <TopoDS.hxx>
94 #include <TopoDS_Shape.hxx>
95 #include <TopoDS_Solid.hxx>
96 #include <TopoDS_Face.hxx>
97 #include <TopoDS_Edge.hxx>
98 #include <TopoDS_Vertex.hxx>
99 #include <TopoDS_Compound.hxx>
100 #include <TopoDS_Iterator.hxx>
101 #include <TopTools_Array1OfShape.hxx>
102 #include <TopTools_DataMapIteratorOfDataMapOfShapeListOfShape.hxx>
103 #include <TopTools_IndexedMapOfShape.hxx>
104 #include <TopTools_ListIteratorOfListOfShape.hxx>
105 #include <TopTools_MapOfShape.hxx>
106 #include <TopTools_MapOfOrientedShape.hxx>
107
108 #include <Geom_Surface.hxx>
109 #include <Geom_Plane.hxx>
110 #include <Geom_SphericalSurface.hxx>
111 #include <Geom_CylindricalSurface.hxx>
112 #include <GeomAdaptor_Surface.hxx>
113
114 #include <GeomLib_Tool.hxx>
115 #include <Geom2d_Curve.hxx>
116
117 #include <Bnd_Box.hxx>
118 #include <GProp_GProps.hxx>
119 #include <TColStd_Array1OfReal.hxx>
120 #include <TColStd_HArray1OfInteger.hxx>
121 #include <TColStd_ListIteratorOfListOfInteger.hxx>
122 #include <TColStd_ListOfInteger.hxx>
123 #include <gp_Cylinder.hxx>
124 #include <gp_Lin.hxx>
125 #include <gp_Pnt.hxx>
126
127 #include <vector>
128 #include <algorithm>
129 #include <functional>
130
131 #include <Standard_NullObject.hxx>
132 #include <Standard_Failure.hxx>
133 #include <Standard_ErrorHandler.hxx> // CAREFUL ! position of this file is critic : see Lucien PIGNOLONI / OCC
134
135 // Includes added for GetInPlace algorithm improvement
136
137 #include <GEOMImpl_MeasureDriver.hxx>
138 #include <GEOMImpl_IMeasure.hxx>
139 #include <BRepBuilderAPI_MakeVertex.hxx>
140
141 #include <BRepClass_FaceClassifier.hxx>
142 #include <BRepClass3d_SolidClassifier.hxx>
143 #include <Precision.hxx>
144
145 #define STD_SORT_ALGO 1
146
147 //=============================================================================
148 /*!
149  *   constructor:
150  */
151 //=============================================================================
152 GEOMImpl_IShapesOperations::GEOMImpl_IShapesOperations (GEOM_Engine* theEngine, int theDocID)
153 : GEOM_IOperations(theEngine, theDocID)
154 {
155   MESSAGE("GEOMImpl_IShapesOperations::GEOMImpl_IShapesOperations");
156 }
157
158 //=============================================================================
159 /*!
160  *  destructor
161  */
162 //=============================================================================
163 GEOMImpl_IShapesOperations::~GEOMImpl_IShapesOperations()
164 {
165   MESSAGE("GEOMImpl_IShapesOperations::~GEOMImpl_IShapesOperations");
166 }
167
168 //=============================================================================
169 /*!
170  *  MakeEdge
171  */
172 //=============================================================================
173 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeEdge
174                      (Handle(GEOM_Object) thePnt1, Handle(GEOM_Object) thePnt2)
175 {
176   SetErrorCode(KO);
177
178   if (thePnt1.IsNull() || thePnt2.IsNull()) return NULL;
179
180   //Add a new Edge object
181   Handle(GEOM_Object) anEdge = GetEngine()->AddObject(GetDocID(), GEOM_EDGE);
182
183   //Add a new Vector function
184   Handle(GEOM_Function) aFunction =
185     anEdge->AddFunction(GEOMImpl_VectorDriver::GetID(), VECTOR_TWO_PNT);
186
187   //Check if the function is set correctly
188   if (aFunction->GetDriverGUID() != GEOMImpl_VectorDriver::GetID()) return NULL;
189
190   GEOMImpl_IVector aPI (aFunction);
191
192   Handle(GEOM_Function) aRef1 = thePnt1->GetLastFunction();
193   Handle(GEOM_Function) aRef2 = thePnt2->GetLastFunction();
194   if (aRef1.IsNull() || aRef2.IsNull()) return NULL;
195
196   aPI.SetPoint1(aRef1);
197   aPI.SetPoint2(aRef2);
198
199   //Compute the Edge value
200   try {
201 #if OCC_VERSION_LARGE > 0x06010000
202     OCC_CATCH_SIGNALS;
203 #endif
204     if (!GetSolver()->ComputeFunction(aFunction)) {
205       SetErrorCode("Vector driver failed");
206       return NULL;
207     }
208   }
209   catch (Standard_Failure) {
210     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
211     SetErrorCode(aFail->GetMessageString());
212     return NULL;
213   }
214
215   //Make a Python command
216   GEOM::TPythonDump(aFunction) << anEdge << " = geompy.MakeEdge("
217                                << thePnt1 << ", " << thePnt2 << ")";
218
219   SetErrorCode(OK);
220   return anEdge;
221 }
222
223 //=============================================================================
224 /*!
225  *  MakeEdgeOnCurveByLength
226  */
227 //=============================================================================
228 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeEdgeOnCurveByLength
229                      (Handle(GEOM_Object) theRefCurve,
230                       const Standard_Real theLength,
231                       Handle(GEOM_Object) theStartPoint)
232 {
233   SetErrorCode(KO);
234
235   if (theRefCurve.IsNull()) return NULL;
236
237   //Add a new Edge object
238   Handle(GEOM_Object) anEdge = GetEngine()->AddObject(GetDocID(), GEOM_EDGE);
239
240   //Add a new Vector function
241   Handle(GEOM_Function) aFunction =
242     anEdge->AddFunction(GEOMImpl_ShapeDriver::GetID(), EDGE_CURVE_LENGTH);
243
244   //Check if the function is set correctly
245   if (aFunction->GetDriverGUID() != GEOMImpl_ShapeDriver::GetID()) return NULL;
246
247   GEOMImpl_IVector aPI (aFunction);
248
249   Handle(GEOM_Function) aRef1 = theRefCurve->GetLastFunction();
250   if (aRef1.IsNull()) return NULL;
251   aPI.SetPoint1(aRef1);
252
253   if (!theStartPoint.IsNull()) {
254     Handle(GEOM_Function) aRef2 = theStartPoint->GetLastFunction();
255     aPI.SetPoint2(aRef2);
256   }
257
258   aPI.SetParameter(theLength);
259
260   //Compute the Edge value
261   try {
262 #if OCC_VERSION_LARGE > 0x06010000
263     OCC_CATCH_SIGNALS;
264 #endif
265     if (!GetSolver()->ComputeFunction(aFunction)) {
266       SetErrorCode("Vector driver failed");
267       return NULL;
268     }
269   }
270   catch (Standard_Failure) {
271     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
272     SetErrorCode(aFail->GetMessageString());
273     return NULL;
274   }
275
276   //Make a Python command
277   GEOM::TPythonDump(aFunction) << anEdge << " = geompy.MakeEdgeOnCurveByLength("
278                                << theRefCurve << ", " << theLength << ", " << theStartPoint << ")";
279
280   SetErrorCode(OK);
281   return anEdge;
282 }
283
284 //=============================================================================
285 /*!
286  *  MakeEdgeWire
287  */
288 //=============================================================================
289 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeEdgeWire
290                     (Handle(GEOM_Object) theWire,
291                      const Standard_Real theLinearTolerance,
292                      const Standard_Real theAngularTolerance)
293 {
294   SetErrorCode(KO);
295
296   if (theWire.IsNull()) return NULL;
297
298   //Add a new Edge object
299   Handle(GEOM_Object) anEdge = GetEngine()->AddObject(GetDocID(), GEOM_EDGE);
300
301   //Add a new Vector function
302   Handle(GEOM_Function) aFunction =
303     anEdge->AddFunction(GEOMImpl_ShapeDriver::GetID(), EDGE_WIRE);
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) aWire = theWire->GetLastFunction();
311
312   if (aWire.IsNull()) return NULL;
313
314   aCI.SetBase(aWire);
315   aCI.SetTolerance(theLinearTolerance);
316   aCI.SetAngularTolerance(theAngularTolerance);
317
318   //Compute the Edge value
319   try {
320 #if OCC_VERSION_LARGE > 0x06010000
321     OCC_CATCH_SIGNALS;
322 #endif
323     if (!GetSolver()->ComputeFunction(aFunction)) {
324       SetErrorCode("Shape driver failed");
325       return NULL;
326     }
327   }
328   catch (Standard_Failure) {
329     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
330     SetErrorCode(aFail->GetMessageString());
331     return NULL;
332   }
333
334   const double DEF_LIN_TOL = Precision::Confusion();
335   const double DEF_ANG_TOL = Precision::Angular();
336   //Make a Python command
337   if ( theAngularTolerance == DEF_ANG_TOL ) {
338     if ( theLinearTolerance == DEF_LIN_TOL )
339       GEOM::TPythonDump(aFunction) << anEdge  << " = geompy.MakeEdgeWire("
340                                    << theWire << ")";
341     else
342       GEOM::TPythonDump(aFunction) << anEdge  << " = geompy.MakeEdgeWire("
343                                    << theWire << ", " << theLinearTolerance << ")";
344   }
345   else {
346     GEOM::TPythonDump(aFunction) << anEdge  << " = geompy.MakeEdgeWire("
347                                  << theWire << ", " << theLinearTolerance << ", "
348                                  << theAngularTolerance << ")";
349   }
350
351   SetErrorCode(OK);
352   return anEdge;
353 }
354
355 //=============================================================================
356 /*!
357  *  MakeWire
358  */
359 //=============================================================================
360 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeWire
361                              (std::list<Handle(GEOM_Object)> theShapes,
362                               const Standard_Real            theTolerance)
363 {
364   SetErrorCode(KO);
365
366   //Add a new object
367   Handle(GEOM_Object) aWire = GetEngine()->AddObject(GetDocID(), GEOM_WIRE);
368
369   //Add a new function
370   Handle(GEOM_Function) aFunction =
371     aWire->AddFunction(GEOMImpl_ShapeDriver::GetID(), WIRE_EDGES);
372   if (aFunction.IsNull()) return NULL;
373
374   //Check if the function is set correctly
375   if (aFunction->GetDriverGUID() != GEOMImpl_ShapeDriver::GetID()) return NULL;
376
377   GEOMImpl_IShapes aCI (aFunction);
378   aCI.SetTolerance(theTolerance);
379
380   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
381
382   // Shapes
383   std::list<Handle(GEOM_Object)>::iterator it = theShapes.begin();
384   for (; it != theShapes.end(); it++) {
385     Handle(GEOM_Function) aRefSh = (*it)->GetLastFunction();
386     if (aRefSh.IsNull()) {
387       SetErrorCode("NULL argument shape for the shape construction");
388       return NULL;
389     }
390     aShapesSeq->Append(aRefSh);
391   }
392   aCI.SetShapes(aShapesSeq);
393
394   //Compute the shape
395   try {
396 #if OCC_VERSION_LARGE > 0x06010000
397     OCC_CATCH_SIGNALS;
398 #endif
399     if (!GetSolver()->ComputeFunction(aFunction)) {
400       SetErrorCode("Shape driver failed");
401       return NULL;
402     }
403   }
404   catch (Standard_Failure) {
405     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
406     SetErrorCode(aFail->GetMessageString());
407     return NULL;
408   }
409
410   //Make a Python command
411   GEOM::TPythonDump pd (aFunction);
412   pd << aWire << " = geompy.MakeWire([";
413
414   // Shapes
415   it = theShapes.begin();
416   if (it != theShapes.end()) {
417     pd << (*it++);
418     while (it != theShapes.end()) {
419       pd << ", " << (*it++);
420     }
421   }
422   pd << "], " << theTolerance << ")";
423
424   SetErrorCode(OK);
425   return aWire;
426 }
427
428 //=============================================================================
429 /*!
430  *  MakeFace
431  */
432 //=============================================================================
433 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeFace (Handle(GEOM_Object) theWire,
434                                                           const bool isPlanarWanted)
435 {
436   SetErrorCode(KO);
437
438   if (theWire.IsNull()) return NULL;
439
440   //Add a new Face object
441   Handle(GEOM_Object) aFace = GetEngine()->AddObject(GetDocID(), GEOM_FACE);
442
443   //Add a new Shape function for creation of a face from a wire
444   Handle(GEOM_Function) aFunction =
445     aFace->AddFunction(GEOMImpl_ShapeDriver::GetID(), FACE_WIRE);
446   if (aFunction.IsNull()) return NULL;
447
448   //Check if the function is set correctly
449   if (aFunction->GetDriverGUID() != GEOMImpl_ShapeDriver::GetID()) return NULL;
450
451   GEOMImpl_IShapes aCI (aFunction);
452
453   Handle(GEOM_Function) aRefWire = theWire->GetLastFunction();
454
455   if (aRefWire.IsNull()) return NULL;
456
457   aCI.SetBase(aRefWire);
458   aCI.SetIsPlanar(isPlanarWanted);
459
460   //Compute the Face value
461   Standard_Boolean isWarning = Standard_False;
462   try {
463 #if OCC_VERSION_LARGE > 0x06010000
464     OCC_CATCH_SIGNALS;
465 #endif
466     if (!GetSolver()->ComputeFunction(aFunction)) {
467       SetErrorCode("Shape driver failed to compute a face");
468       return NULL;
469     }
470   }
471   catch (Standard_Failure) {
472     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
473     SetErrorCode(aFail->GetMessageString());
474     // to provide warning
475     if (!aFunction->GetValue().IsNull()) {
476       isWarning = Standard_True;
477     } else {
478       return NULL;
479     }
480   }
481
482   //Make a Python command
483   GEOM::TPythonDump(aFunction) << aFace << " = geompy.MakeFace("
484     << theWire << ", " << (int)isPlanarWanted << ")";
485
486   // to provide warning
487   if (!isWarning) SetErrorCode(OK);
488   return aFace;
489 }
490
491 //=============================================================================
492 /*!
493  *  MakeFaceWires
494  */
495 //=============================================================================
496 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeFaceWires
497                              (std::list<Handle(GEOM_Object)> theShapes,
498                               const bool isPlanarWanted)
499 {
500   SetErrorCode(KO);
501
502   //Add a new object
503   Handle(GEOM_Object) aShape = GetEngine()->AddObject(GetDocID(), GEOM_FACE);
504
505   //Add a new function
506   Handle(GEOM_Function) aFunction =
507     aShape->AddFunction(GEOMImpl_ShapeDriver::GetID(), FACE_WIRES);
508   if (aFunction.IsNull()) return NULL;
509
510   //Check if the function is set correctly
511   if (aFunction->GetDriverGUID() != GEOMImpl_ShapeDriver::GetID()) return NULL;
512
513   GEOMImpl_IShapes aCI (aFunction);
514
515   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
516
517   // Shapes
518   std::list<Handle(GEOM_Object)>::iterator it = theShapes.begin();
519   for (; it != theShapes.end(); it++) {
520     Handle(GEOM_Function) aRefSh = (*it)->GetLastFunction();
521     if (aRefSh.IsNull()) {
522       SetErrorCode("NULL argument shape for the face construction");
523       return NULL;
524     }
525     aShapesSeq->Append(aRefSh);
526   }
527   aCI.SetShapes(aShapesSeq);
528
529   aCI.SetIsPlanar(isPlanarWanted);
530
531   //Compute the shape
532   Standard_Boolean isWarning = Standard_False;
533   try {
534 #if OCC_VERSION_LARGE > 0x06010000
535     OCC_CATCH_SIGNALS;
536 #endif
537     if (!GetSolver()->ComputeFunction(aFunction)) {
538       SetErrorCode("Shape driver failed");
539       return NULL;
540     }
541   }
542   catch (Standard_Failure) {
543     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
544     SetErrorCode(aFail->GetMessageString());
545     // to provide warning
546     if (!aFunction->GetValue().IsNull()) {
547       isWarning = Standard_True;
548     } else {
549       return NULL;
550     }
551   }
552
553   //Make a Python command
554   GEOM::TPythonDump pd (aFunction);
555   pd << aShape << " = geompy.MakeFaceWires([";
556
557   // Shapes
558   it = theShapes.begin();
559   if (it != theShapes.end()) {
560     pd << (*it++);
561     while (it != theShapes.end()) {
562       pd << ", " << (*it++);
563     }
564   }
565   pd << "], " << (int)isPlanarWanted << ")";
566
567   // to provide warning
568   if (!isWarning) SetErrorCode(OK);
569   return aShape;
570 }
571
572 //=============================================================================
573 /*!
574  *  MakeShell
575  */
576 //=============================================================================
577 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeShell
578                              (std::list<Handle(GEOM_Object)> theShapes)
579 {
580   return MakeShape(theShapes, GEOM_SHELL, SHELL_FACES, "MakeShell");
581 }
582
583 //=============================================================================
584 /*!
585  *  MakeSolidShells
586  */
587 //=============================================================================
588 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeSolidShells
589                              (std::list<Handle(GEOM_Object)> theShapes)
590 {
591   return MakeShape(theShapes, GEOM_SOLID, SOLID_SHELLS, "MakeSolid");
592 }
593
594 //=============================================================================
595 /*!
596  *  MakeCompound
597  */
598 //=============================================================================
599 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeCompound
600                              (std::list<Handle(GEOM_Object)> theShapes)
601 {
602   return MakeShape(theShapes, GEOM_COMPOUND, COMPOUND_SHAPES, "MakeCompound");
603 }
604
605 //=============================================================================
606 /*!
607  *  MakeShape
608  */
609 //=============================================================================
610 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeShape
611                              (std::list<Handle(GEOM_Object)> theShapes,
612                               const Standard_Integer         theObjectType,
613                               const Standard_Integer         theFunctionType,
614                               const TCollection_AsciiString& theMethodName)
615 {
616   SetErrorCode(KO);
617
618   //Add a new object
619   Handle(GEOM_Object) aShape = GetEngine()->AddObject(GetDocID(), theObjectType);
620
621   //Add a new function
622   Handle(GEOM_Function) aFunction =
623     aShape->AddFunction(GEOMImpl_ShapeDriver::GetID(), theFunctionType);
624   if (aFunction.IsNull()) return NULL;
625
626   //Check if the function is set correctly
627   if (aFunction->GetDriverGUID() != GEOMImpl_ShapeDriver::GetID()) return NULL;
628
629   GEOMImpl_IShapes aCI (aFunction);
630
631   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
632
633   // Shapes
634   std::list<Handle(GEOM_Object)>::iterator it = theShapes.begin();
635   for (; it != theShapes.end(); it++) {
636     Handle(GEOM_Function) aRefSh = (*it)->GetLastFunction();
637     if (aRefSh.IsNull()) {
638       SetErrorCode("NULL argument shape for the shape construction");
639       return NULL;
640     }
641     aShapesSeq->Append(aRefSh);
642   }
643   aCI.SetShapes(aShapesSeq);
644
645   //Compute the shape
646   try {
647 #if OCC_VERSION_LARGE > 0x06010000
648     OCC_CATCH_SIGNALS;
649 #endif
650     if (!GetSolver()->ComputeFunction(aFunction)) {
651       SetErrorCode("Shape driver failed");
652       return NULL;
653     }
654   }
655   catch (Standard_Failure) {
656     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
657     SetErrorCode(aFail->GetMessageString());
658     return NULL;
659   }
660
661   //Make a Python command
662   GEOM::TPythonDump pd (aFunction);
663   pd << aShape << " = geompy." << theMethodName.ToCString() << "([";
664
665   // Shapes
666   it = theShapes.begin();
667   if (it != theShapes.end()) {
668     pd << (*it++);
669     while (it != theShapes.end()) {
670       pd << ", " << (*it++);
671     }
672   }
673   pd << "])";
674
675   SetErrorCode(OK);
676   return aShape;
677 }
678
679 //=============================================================================
680 /*!
681  *  MakeGlueFaces
682  */
683 //=============================================================================
684 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeGlueFaces
685                                                 (Handle(GEOM_Object) theShape,
686                                                  const Standard_Real theTolerance,
687                                                  const Standard_Boolean doKeepNonSolids)
688 {
689   SetErrorCode(KO);
690
691   if (theShape.IsNull()) return NULL;
692
693   //Add a new Glued object
694   Handle(GEOM_Object) aGlued = GetEngine()->AddObject(GetDocID(), GEOM_GLUED);
695
696   //Add a new Glue function
697   Handle(GEOM_Function) aFunction;
698   aFunction = aGlued->AddFunction(GEOMImpl_GlueDriver::GetID(), GLUE_FACES);
699   if (aFunction.IsNull()) return NULL;
700
701   //Check if the function is set correctly
702   if (aFunction->GetDriverGUID() != GEOMImpl_GlueDriver::GetID()) return NULL;
703
704   GEOMImpl_IGlue aCI (aFunction);
705
706   Handle(GEOM_Function) aRefShape = theShape->GetLastFunction();
707   if (aRefShape.IsNull()) return NULL;
708
709   aCI.SetBase(aRefShape);
710   aCI.SetTolerance(theTolerance);
711   aCI.SetKeepNonSolids(doKeepNonSolids);
712
713   //Compute the sub-shape value
714   Standard_Boolean isWarning = Standard_False;
715   try {
716 #if OCC_VERSION_LARGE > 0x06010000
717     OCC_CATCH_SIGNALS;
718 #endif
719     if (!GetSolver()->ComputeFunction(aFunction)) {
720       SetErrorCode("Shape driver failed to glue faces");
721       return NULL;
722     }
723   }
724   catch (Standard_Failure) {
725     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
726     SetErrorCode(aFail->GetMessageString());
727     // to provide warning
728     if (!aFunction->GetValue().IsNull()) {
729       isWarning = Standard_True;
730     } else {
731       return NULL;
732     }
733   }
734
735   //Make a Python command
736   GEOM::TPythonDump(aFunction) << aGlued << " = geompy.MakeGlueFaces("
737     << theShape << ", " << theTolerance << ")";
738
739   // to provide warning
740   if (!isWarning) SetErrorCode(OK);
741   return aGlued;
742 }
743
744 //=============================================================================
745 /*!
746  *  GetGlueFaces
747  */
748 //=============================================================================
749 /*
750 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetGlueFaces
751                                                 (Handle(GEOM_Object) theShape,
752                                                  const Standard_Real theTolerance)
753 {
754   SetErrorCode(KO);
755
756   if (theShape.IsNull()) return NULL;
757   TopoDS_Shape aShape = theShape->GetValue();
758   if (aShape.IsNull()) return NULL;
759
760   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
761
762   Standard_Integer iErr;
763   TopoDS_Shape aS;
764   GEOMAlgo_Gluer1 aGluer;
765   GEOMAlgo_ListIteratorOfListOfCoupleOfShapes aItCS;
766   GEOMAlgo_CoupleOfShapes aCS;
767   GEOMAlgo_ListOfCoupleOfShapes aLCS;
768
769   //aGluer = new GEOMAlgo_Gluer1;
770   aGluer.SetShape(aShape);
771   aGluer.SetTolerance(theTolerance);
772   aGluer.Perform();
773   iErr = aGluer.ErrorStatus();
774   if (iErr) return NULL;
775
776   TopTools_ListOfShape listShape;
777   const GEOMAlgo_ListOfCoupleOfShapes& aLCSG = aGluer.GluedFaces();
778   // Access to faces
779   aItCS.Initialize(aLCSG);
780   for (; aItCS.More(); aItCS.Next()) {
781     const GEOMAlgo_CoupleOfShapes& aCSG = aItCS.Value();
782     listShape.Append(aCSG.Shape1());
783   }
784
785   TopTools_ListIteratorOfListOfShape itSub (listShape);
786   TCollection_AsciiString anAsciiList, anEntry;
787   TopTools_IndexedMapOfShape anIndices;
788   TopExp::MapShapes(aShape, anIndices);
789   Handle(TColStd_HArray1OfInteger) anArray;
790   Handle(GEOM_Object) anObj;
791   for (int index = 1; itSub.More(); itSub.Next(), ++index) {
792     TopoDS_Shape aValue = itSub.Value();
793     anArray = new TColStd_HArray1OfInteger(1,1);
794     anArray->SetValue(1, anIndices.FindIndex(aValue));
795     anObj = GetEngine()->AddSubShape(theShape, anArray);
796     if (!anObj.IsNull()) {
797       aSeq->Append(anObj);
798
799       // for python command
800       TDF_Tool::Entry(anObj->GetEntry(), anEntry);
801       anAsciiList += anEntry;
802       anAsciiList += ",";
803     }
804   }
805
806   //Make a Python command
807   if( anAsciiList.Length() > 0 ) {
808     anAsciiList.Trunc(anAsciiList.Length() - 1);
809     Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
810     GEOM::TPythonDump pd (aFunction, true);
811     pd << "[" << anAsciiList.ToCString();
812     pd << "] = geompy.GetGlueFaces(" << theShape << ", " << theTolerance << ")";
813   }
814
815   SetErrorCode(OK);
816
817   return aSeq;
818 }
819 */
820
821 //=============================================================================
822 /*!
823  *  MakeGlueFacesByList
824  */
825 //=============================================================================
826 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeGlueFacesByList
827                                                 (Handle(GEOM_Object) theShape,
828                                                  const Standard_Real theTolerance,
829                                                  std::list<Handle(GEOM_Object)> theFaces,
830                                                  const Standard_Boolean doKeepNonSolids,
831                                                  const Standard_Boolean doGlueAllEdges)
832 {
833   SetErrorCode(KO);
834
835   if (theShape.IsNull()) return NULL;
836
837   //Add a new Glued object
838   Handle(GEOM_Object) aGlued = GetEngine()->AddObject(GetDocID(), GEOM_GLUED);
839
840   //Add a new Glue function
841   Handle(GEOM_Function) aFunction;
842   aFunction = aGlued->AddFunction(GEOMImpl_GlueDriver::GetID(), GLUE_FACES_BY_LIST);
843   if (aFunction.IsNull()) return NULL;
844
845   //Check if the function is set correctly
846   if (aFunction->GetDriverGUID() != GEOMImpl_GlueDriver::GetID()) return NULL;
847
848   GEOMImpl_IGlue aCI (aFunction);
849
850   Handle(GEOM_Function) aRefShape = theShape->GetLastFunction();
851   if (aRefShape.IsNull()) return NULL;
852
853   aCI.SetBase(aRefShape);
854   aCI.SetTolerance(theTolerance);
855   aCI.SetKeepNonSolids(doKeepNonSolids);
856   aCI.SetGlueAllEdges(doGlueAllEdges);
857
858   Handle(TColStd_HSequenceOfTransient) aFaces = new TColStd_HSequenceOfTransient;
859   std::list<Handle(GEOM_Object)>::iterator it = theFaces.begin();
860   for (; it != theFaces.end(); it++) {
861     Handle(GEOM_Function) aRefSh = (*it)->GetLastFunction();
862     if (aRefSh.IsNull()) {
863       SetErrorCode("NULL argument shape for the shape construction");
864       return NULL;
865     }
866     aFaces->Append(aRefSh);
867   }
868   aCI.SetFaces(aFaces);
869
870   //Compute the sub-shape value
871   Standard_Boolean isWarning = Standard_False;
872   try {
873 #if OCC_VERSION_LARGE > 0x06010000
874     OCC_CATCH_SIGNALS;
875 #endif
876     if (!GetSolver()->ComputeFunction(aFunction)) {
877       SetErrorCode("Shape driver failed to glue faces");
878       return NULL;
879     }
880   }
881   catch (Standard_Failure) {
882     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
883     SetErrorCode(aFail->GetMessageString());
884     // to provide warning
885     if (!aFunction->GetValue().IsNull()) {
886       isWarning = Standard_True;
887     } else {
888       return NULL;
889     }
890   }
891
892   //Make a Python command
893
894   GEOM::TPythonDump pd(aFunction);
895   pd << aGlued << " = geompy.MakeGlueFacesByList("
896      << theShape << ", " << theTolerance << ", [";
897   // Faces
898   it = theFaces.begin();
899   if (it != theFaces.end()) {
900     pd << (*it++);
901     while (it != theFaces.end()) {
902       pd << ", " << (*it++);
903     }
904   }
905   pd << "], " << (bool)doKeepNonSolids << ", " << (bool)doGlueAllEdges << ")";
906
907   // to provide warning
908   if (!isWarning) SetErrorCode(OK);
909   return aGlued;
910 }
911
912 //=============================================================================
913 /*!
914  *  MakeGlueEdges
915  */
916 //=============================================================================
917 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeGlueEdges
918                                                 (Handle(GEOM_Object) theShape,
919                                                  const Standard_Real theTolerance)
920 {
921   SetErrorCode(KO);
922
923   if (theShape.IsNull()) return NULL;
924
925   //Add a new Glued object
926   Handle(GEOM_Object) aGlued = GetEngine()->AddObject(GetDocID(), GEOM_GLUED);
927
928   //Add a new Glue function
929   Handle(GEOM_Function) aFunction;
930   aFunction = aGlued->AddFunction(GEOMImpl_GlueDriver::GetID(), GLUE_EDGES);
931   if (aFunction.IsNull()) return NULL;
932
933   //Check if the function is set correctly
934   if (aFunction->GetDriverGUID() != GEOMImpl_GlueDriver::GetID()) return NULL;
935
936   GEOMImpl_IGlue aCI (aFunction);
937
938   Handle(GEOM_Function) aRefShape = theShape->GetLastFunction();
939   if (aRefShape.IsNull()) return NULL;
940
941   aCI.SetBase(aRefShape);
942   aCI.SetTolerance(theTolerance);
943   aCI.SetKeepNonSolids(true);
944
945   //Compute the sub-shape value
946   Standard_Boolean isWarning = Standard_False;
947   try {
948 #if OCC_VERSION_LARGE > 0x06010000
949     OCC_CATCH_SIGNALS;
950 #endif
951     if (!GetSolver()->ComputeFunction(aFunction)) {
952       SetErrorCode("Shape driver failed to glue edges");
953       return NULL;
954     }
955   }
956   catch (Standard_Failure) {
957     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
958     SetErrorCode(aFail->GetMessageString());
959     // to provide warning
960     if (!aFunction->GetValue().IsNull()) {
961       isWarning = Standard_True;
962     } else {
963       return NULL;
964     }
965   }
966
967   //Make a Python command
968   GEOM::TPythonDump(aFunction) << aGlued << " = geompy.MakeGlueEdges("
969     << theShape << ", " << theTolerance << ")";
970
971   // to provide warning
972   if (!isWarning) SetErrorCode(OK);
973   return aGlued;
974 }
975
976 //=============================================================================
977 /*!
978  *  GetGlueShapes
979  */
980 //=============================================================================
981 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetGlueShapes
982                                                 (Handle(GEOM_Object)    theShape,
983                                                  const Standard_Real    theTolerance,
984                                                  const TopAbs_ShapeEnum theType)
985 {
986   SetErrorCode(KO);
987
988   if (theShape.IsNull()) return NULL;
989   TopoDS_Shape aShape = theShape->GetValue();
990   if (aShape.IsNull()) return NULL;
991
992   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
993
994   GEOMAlgo_GlueDetector aGluer;
995   aGluer.SetArgument(aShape);
996   aGluer.SetTolerance(theTolerance);
997   aGluer.Perform();
998   Standard_Integer iErr = aGluer.ErrorStatus();
999   if (iErr) return NULL;
1000
1001   TCollection_AsciiString anAsciiList, anEntry;
1002   TopTools_IndexedMapOfShape anIndices;
1003   TopExp::MapShapes(aShape, anIndices);
1004   Handle(TColStd_HArray1OfInteger) anArray;
1005   Handle(GEOM_Object) anObj;
1006
1007   TopTools_ListOfShape listOnePerSet;
1008
1009   const TopTools_DataMapOfShapeListOfShape& aImages = aGluer.Images();
1010   TopTools_DataMapIteratorOfDataMapOfShapeListOfShape aItDMSLS (aImages);
1011   for (int index = 1; aItDMSLS.More(); aItDMSLS.Next(), ++index) {
1012     // some key shape
1013     //const TopoDS_Shape& aSkey = aItDMSLS.Key();  
1014
1015     // list of shapes of the argument that can be glued
1016     const TopTools_ListOfShape& aLSD = aItDMSLS.Value(); 
1017
1018     //listShape.Append(aLSD.First());
1019     TopoDS_Shape aValue = aLSD.First();
1020
1021     if (aValue.ShapeType() == theType) {
1022       listOnePerSet.Append(aValue);
1023     }
1024   }
1025
1026   // for stable order of returned entities
1027   GEOMImpl_IShapesOperations::SortShapes(listOnePerSet, Standard_False);
1028
1029   TopTools_ListIteratorOfListOfShape aListIt (listOnePerSet);
1030   for (; aListIt.More(); aListIt.Next()) {
1031     TopoDS_Shape aValue = aListIt.Value();
1032     anArray = new TColStd_HArray1OfInteger(1,1);
1033     anArray->SetValue(1, anIndices.FindIndex(aValue));
1034     anObj = GetEngine()->AddSubShape(theShape, anArray);
1035     if (!anObj.IsNull()) {
1036       aSeq->Append(anObj);
1037
1038       // for python command
1039       TDF_Tool::Entry(anObj->GetEntry(), anEntry);
1040       anAsciiList += anEntry;
1041       anAsciiList += ",";
1042     }
1043   }
1044
1045   // Make a Python command
1046   if (anAsciiList.Length() > 0) {
1047     anAsciiList.Trunc(anAsciiList.Length() - 1);
1048     Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
1049     GEOM::TPythonDump pd (aFunction, /*append=*/true);
1050     pd << "[" << anAsciiList.ToCString();
1051     if (theType == TopAbs_FACE)
1052       pd << "] = geompy.GetGlueFaces(" << theShape << ", " << theTolerance << ")";
1053     else if (theType == TopAbs_EDGE)
1054       pd << "] = geompy.GetGlueEdges(" << theShape << ", " << theTolerance << ")";
1055   }
1056
1057   SetErrorCode(OK);
1058
1059   return aSeq;
1060 }
1061
1062 //=============================================================================
1063 /*!
1064  *  MakeGlueEdgesByList
1065  */
1066 //=============================================================================
1067 Handle(GEOM_Object) GEOMImpl_IShapesOperations::MakeGlueEdgesByList
1068                                                 (Handle(GEOM_Object) theShape,
1069                                                  const Standard_Real theTolerance,
1070                                                  std::list<Handle(GEOM_Object)> theEdges)
1071 {
1072   SetErrorCode(KO);
1073
1074   if (theShape.IsNull()) return NULL;
1075
1076   //Add a new Glued object
1077   Handle(GEOM_Object) aGlued = GetEngine()->AddObject(GetDocID(), GEOM_GLUED);
1078
1079   //Add a new Glue function
1080   Handle(GEOM_Function) aFunction;
1081   aFunction = aGlued->AddFunction(GEOMImpl_GlueDriver::GetID(), GLUE_EDGES_BY_LIST);
1082   if (aFunction.IsNull()) return NULL;
1083
1084   //Check if the function is set correctly
1085   if (aFunction->GetDriverGUID() != GEOMImpl_GlueDriver::GetID()) return NULL;
1086
1087   GEOMImpl_IGlue aCI (aFunction);
1088
1089   Handle(GEOM_Function) aRefShape = theShape->GetLastFunction();
1090   if (aRefShape.IsNull()) return NULL;
1091
1092   aCI.SetBase(aRefShape);
1093   aCI.SetTolerance(theTolerance);
1094   aCI.SetKeepNonSolids(true);
1095
1096   Handle(TColStd_HSequenceOfTransient) anEdges = new TColStd_HSequenceOfTransient;
1097   std::list<Handle(GEOM_Object)>::iterator it = theEdges.begin();
1098   for (; it != theEdges.end(); it++) {
1099     Handle(GEOM_Function) aRefSh = (*it)->GetLastFunction();
1100     if (aRefSh.IsNull()) {
1101       SetErrorCode("NULL argument shape for the shape construction");
1102       return NULL;
1103     }
1104     anEdges->Append(aRefSh);
1105   }
1106   aCI.SetFaces(anEdges);
1107
1108   //Compute the sub-shape value
1109   Standard_Boolean isWarning = Standard_False;
1110   try {
1111 #if OCC_VERSION_LARGE > 0x06010000
1112     OCC_CATCH_SIGNALS;
1113 #endif
1114     if (!GetSolver()->ComputeFunction(aFunction)) {
1115       SetErrorCode("Shape driver failed to glue edges");
1116       return NULL;
1117     }
1118   }
1119   catch (Standard_Failure) {
1120     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1121     SetErrorCode(aFail->GetMessageString());
1122     // to provide warning
1123     if (!aFunction->GetValue().IsNull()) {
1124       isWarning = Standard_True;
1125     } else {
1126       return NULL;
1127     }
1128   }
1129
1130   //Make a Python command
1131
1132   GEOM::TPythonDump pd (aFunction);
1133   pd << aGlued << " = geompy.MakeGlueEdgesByList("
1134      << theShape << ", " << theTolerance << ", [";
1135   // Edges
1136   it = theEdges.begin();
1137   if (it != theEdges.end()) {
1138     pd << (*it++);
1139     while (it != theEdges.end()) {
1140       pd << ", " << (*it++);
1141     }
1142   }
1143   pd << "])";
1144
1145   // to provide warning
1146   if (!isWarning) SetErrorCode(OK);
1147   return aGlued;
1148 }
1149
1150 //=============================================================================
1151 /*!
1152  *  GetExistingSubObjects
1153  */
1154 //=============================================================================
1155 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetExistingSubObjects
1156                                           (Handle(GEOM_Object)    theShape,
1157                                            const Standard_Boolean theGroupsOnly)
1158 {
1159   SetErrorCode(KO);
1160
1161   if (theShape.IsNull()) return NULL;
1162
1163   Handle(GEOM_Function) aMainShape = theShape->GetLastFunction();
1164   if (aMainShape.IsNull()) return NULL;
1165
1166   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
1167   SetErrorCode(NOT_FOUND_ANY);
1168
1169   if (!aMainShape->HasSubShapeReferences()) return aSeq;
1170   const TDataStd_ListOfExtendedString& aListEntries = aMainShape->GetSubShapeReferences();
1171   if (aListEntries.IsEmpty()) return aSeq;
1172
1173   SetErrorCode(KO);
1174
1175   TCollection_AsciiString anAsciiList;
1176
1177   TDataStd_ListIteratorOfListOfExtendedString anIt (aListEntries);
1178   for (; anIt.More(); anIt.Next()) {
1179     TCollection_ExtendedString anEntry = anIt.Value();
1180     Standard_Integer aStrLen = anEntry.LengthOfCString();
1181     char* anEntryStr = new char[aStrLen];
1182     anEntry.ToUTF8CString(anEntryStr);
1183     Handle(GEOM_Object) anObj = GetEngine()->GetObject(GetDocID(), anEntryStr, false);
1184     if (!anObj.IsNull()) {
1185       if (!theGroupsOnly || anObj->GetType() == GEOM_GROUP) {
1186         aSeq->Append(anObj);
1187
1188         // for python command
1189         anAsciiList += anEntryStr;
1190         anAsciiList += ",";
1191       }
1192     }
1193     delete [] anEntryStr;
1194   }
1195
1196   if (aSeq->Length() == 0) {
1197     SetErrorCode(NOT_FOUND_ANY);
1198     return aSeq;
1199   }
1200
1201   //Make a Python command
1202   anAsciiList.Trunc(anAsciiList.Length() - 1);
1203
1204   GEOM::TPythonDump pd (aMainShape, /*append=*/true);
1205   pd << "[" << anAsciiList.ToCString();
1206   pd << "] = geompy.GetExistingSubObjects(";
1207   pd << theShape << ", " << (bool)theGroupsOnly << ")";
1208
1209   SetErrorCode(OK);
1210
1211   return aSeq;
1212 }
1213
1214 //=============================================================================
1215 /*!
1216  *  MakeExplode
1217  */
1218 //=============================================================================
1219 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::MakeExplode
1220                                           (Handle(GEOM_Object)    theShape,
1221                                            const Standard_Integer theShapeType,
1222                                            const Standard_Boolean isSorted,
1223                                            const ExplodeType      theExplodeType)
1224 {
1225   SetErrorCode(KO);
1226
1227   if (theShape.IsNull()) return NULL;
1228   TopoDS_Shape aShape = theShape->GetValue();
1229   if (aShape.IsNull()) return NULL;
1230
1231   Handle(GEOM_Function) aMainShape = theShape->GetLastFunction();
1232
1233   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
1234   Handle(GEOM_Object) anObj;
1235   TopTools_MapOfShape mapShape;
1236   TopTools_ListOfShape listShape;
1237
1238   if (aShape.ShapeType() == TopAbs_COMPOUND &&
1239       (TopAbs_ShapeEnum(theShapeType) == TopAbs_SHAPE ||
1240        TopAbs_ShapeEnum(theShapeType) == TopAbs_COMPSOLID ||
1241        TopAbs_ShapeEnum(theShapeType) == TopAbs_COMPOUND))
1242   {
1243     TopoDS_Iterator It (aShape, Standard_True, Standard_True);
1244     for (; It.More(); It.Next()) {
1245       if (mapShape.Add(It.Value())) {
1246         if (TopAbs_ShapeEnum(theShapeType) == TopAbs_SHAPE ||
1247             TopAbs_ShapeEnum(theShapeType) == It.Value().ShapeType()) {
1248           listShape.Append(It.Value());
1249         }
1250       }
1251     }
1252   }
1253   else if (theExplodeType != EXPLODE_NEW_EXCLUDE_MAIN || aShape.ShapeType() != theShapeType) // issue 0021079
1254   {
1255     TopExp_Explorer exp (aShape, TopAbs_ShapeEnum(theShapeType));
1256     for (; exp.More(); exp.Next())
1257       if (mapShape.Add(exp.Current()))
1258         listShape.Append(exp.Current());
1259   }
1260
1261   if (listShape.IsEmpty()) {
1262     //SetErrorCode("The given shape has no sub-shapes of the requested type");
1263     SetErrorCode(NOT_FOUND_ANY); // NPAL18017
1264     return aSeq;
1265   }
1266
1267   if (isSorted) {
1268     bool isOldSorting = false;
1269     if (theExplodeType == EXPLODE_OLD_INCLUDE_MAIN)
1270       isOldSorting = true;
1271     SortShapes(listShape, isOldSorting);
1272   }
1273
1274   TopTools_IndexedMapOfShape anIndices;
1275   TopExp::MapShapes(aShape, anIndices);
1276   Handle(TColStd_HArray1OfInteger) anArray;
1277
1278   TopTools_ListIteratorOfListOfShape itSub (listShape);
1279   TCollection_AsciiString anAsciiList, anEntry;
1280   for (int index = 1; itSub.More(); itSub.Next(), ++index)
1281   {
1282     TopoDS_Shape aValue = itSub.Value();
1283     anArray = new TColStd_HArray1OfInteger(1,1);
1284     anArray->SetValue(1, anIndices.FindIndex(aValue));
1285
1286     //anObj = GetEngine()->AddSubShape(theShape, anArray);
1287     {
1288       anObj = GetEngine()->AddObject(GetDocID(), GEOM_SUBSHAPE);
1289       Handle(GEOM_Function) aFunction = anObj->AddFunction(GEOM_Object::GetSubShapeID(), 1);
1290       if (aFunction.IsNull()) return aSeq;
1291
1292       GEOM_ISubShape aSSI (aFunction);
1293       aSSI.SetMainShape(aMainShape);
1294       aSSI.SetIndices(anArray);
1295
1296       // Set function value directly, as we know it.
1297       // Usage of Solver here would lead to significant loss of time,
1298       // because GEOM_SubShapeDriver will build TopTools_IndexedMapOfShape
1299       // on the main shape for each being calculated sub-shape separately.
1300       aFunction->SetValue(aValue);
1301
1302       // Put this subshape in the list of sub-shapes of theMainShape
1303       aMainShape->AddSubShapeReference(aFunction);
1304     }
1305
1306     if (!anObj.IsNull()) {
1307       aSeq->Append(anObj);
1308
1309       // for python command
1310       TDF_Tool::Entry(anObj->GetEntry(), anEntry);
1311       anAsciiList += anEntry;
1312       anAsciiList += ",";
1313     }
1314   }
1315
1316   //Make a Python command
1317   anAsciiList.Trunc(anAsciiList.Length() - 1);
1318
1319   GEOM::TPythonDump pd (aMainShape, /*append=*/true);
1320   pd << "[" << anAsciiList.ToCString() << "] = geompy.";
1321   switch (theExplodeType) {
1322   case EXPLODE_NEW_EXCLUDE_MAIN:
1323     pd << "ExtractShapes(" << theShape << ", "
1324        << TopAbs_ShapeEnum(theShapeType) << ", " << (isSorted ? "True" : "False") << ")";
1325     break;
1326   case EXPLODE_NEW_INCLUDE_MAIN:
1327     pd << "SubShapeAll" << (isSorted ? "SortedCentres(" : "(")
1328        << theShape << ", " << TopAbs_ShapeEnum(theShapeType) << ")";
1329     break;
1330   case EXPLODE_OLD_INCLUDE_MAIN:
1331     pd << "SubShapeAll" << (isSorted ? "Sorted(" : "(")
1332        << theShape << ", " << TopAbs_ShapeEnum(theShapeType) << ")";
1333     break;
1334   default: ;
1335   }
1336   SetErrorCode(OK);
1337
1338   return aSeq;
1339 }
1340
1341 //=============================================================================
1342 /*!
1343  *  SubShapeAllIDs
1344  */
1345 //=============================================================================
1346 Handle(TColStd_HSequenceOfInteger) GEOMImpl_IShapesOperations::SubShapeAllIDs
1347                                           (Handle(GEOM_Object)    theShape,
1348                                            const Standard_Integer theShapeType,
1349                                            const Standard_Boolean isSorted,
1350                                            const ExplodeType      theExplodeType)
1351 {
1352   SetErrorCode(KO);
1353
1354   if (theShape.IsNull()) return NULL;
1355   TopoDS_Shape aShape = theShape->GetValue();
1356   if (aShape.IsNull()) return NULL;
1357
1358   Handle(TColStd_HSequenceOfInteger) aSeq = new TColStd_HSequenceOfInteger;
1359   TopTools_MapOfShape mapShape;
1360   TopTools_ListOfShape listShape;
1361
1362   if (aShape.ShapeType() == TopAbs_COMPOUND &&
1363       (TopAbs_ShapeEnum(theShapeType) == TopAbs_SHAPE ||
1364        TopAbs_ShapeEnum(theShapeType) == TopAbs_COMPSOLID ||
1365        TopAbs_ShapeEnum(theShapeType) == TopAbs_COMPOUND))
1366   {
1367     TopoDS_Iterator It (aShape, Standard_True, Standard_True);
1368     for (; It.More(); It.Next()) {
1369       if (mapShape.Add(It.Value())) {
1370         if (TopAbs_ShapeEnum(theShapeType) == TopAbs_SHAPE ||
1371             TopAbs_ShapeEnum(theShapeType) == It.Value().ShapeType()) {
1372           listShape.Append(It.Value());
1373         }
1374       }
1375     }
1376   }
1377   else if (theExplodeType != EXPLODE_NEW_EXCLUDE_MAIN || aShape.ShapeType() != theShapeType) // issue 0021079
1378   {
1379     TopExp_Explorer exp (aShape, TopAbs_ShapeEnum(theShapeType));
1380     for (; exp.More(); exp.Next())
1381       if (mapShape.Add(exp.Current()))
1382         listShape.Append(exp.Current());
1383   }
1384
1385   if (listShape.IsEmpty()) {
1386     //SetErrorCode("The given shape has no sub-shapes of the requested type");
1387     SetErrorCode(NOT_FOUND_ANY); // NPAL18017
1388     return aSeq;
1389   }
1390
1391   if (isSorted) {
1392     bool isOldSorting = false;
1393     if (theExplodeType == EXPLODE_OLD_INCLUDE_MAIN)
1394       isOldSorting = true;
1395     SortShapes(listShape, isOldSorting);
1396   }
1397
1398   TopTools_IndexedMapOfShape anIndices;
1399   TopExp::MapShapes(aShape, anIndices);
1400   Handle(TColStd_HArray1OfInteger) anArray;
1401
1402   TopTools_ListIteratorOfListOfShape itSub (listShape);
1403   for (int index = 1; itSub.More(); itSub.Next(), ++index) {
1404     TopoDS_Shape aValue = itSub.Value();
1405     aSeq->Append(anIndices.FindIndex(aValue));
1406   }
1407
1408   Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
1409
1410   //Make a Python command
1411   GEOM::TPythonDump pd (aFunction, /*append=*/true);
1412   pd << "listSubShapeIDs = geompy.SubShapeAll";
1413   switch (theExplodeType) {
1414   case EXPLODE_NEW_EXCLUDE_MAIN:
1415     break;
1416   case EXPLODE_NEW_INCLUDE_MAIN:
1417     pd << (isSorted ? "SortedCentresIDs(" : "IDs(")
1418        << theShape << ", " << TopAbs_ShapeEnum(theShapeType) << ")";
1419     break;
1420   case EXPLODE_OLD_INCLUDE_MAIN:
1421     pd << (isSorted ? "SortedIDs(" : "IDs(")
1422        << theShape << ", " << TopAbs_ShapeEnum(theShapeType) << ")";
1423     break;
1424   default: ;
1425   }
1426
1427   SetErrorCode(OK);
1428   return aSeq;
1429 }
1430
1431 //=============================================================================
1432 /*!
1433  *  GetSubShape
1434  */
1435 //=============================================================================
1436 Handle(GEOM_Object) GEOMImpl_IShapesOperations::GetSubShape
1437                                           (Handle(GEOM_Object)    theMainShape,
1438                                            const Standard_Integer theID)
1439 {
1440   SetErrorCode(KO);
1441
1442   if (theMainShape.IsNull()) return NULL;
1443
1444   Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
1445   anArray->SetValue(1, theID);
1446   Handle(GEOM_Object) anObj = GetEngine()->AddSubShape(theMainShape, anArray,true);
1447   if (anObj.IsNull()) {
1448     SetErrorCode("Can not get a sub-shape with the given ID");
1449     return NULL;
1450   }
1451
1452   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
1453
1454   //Make a Python command
1455   GEOM::TPythonDump(aFunction) << anObj << " = geompy.GetSubShape("
1456                                << theMainShape << ", [" << theID << "])";
1457
1458   SetErrorCode(OK);
1459   return anObj;
1460 }
1461
1462 //=============================================================================
1463 /*!
1464  *  MakeSubShapes
1465  */
1466 //=============================================================================
1467 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::MakeSubShapes
1468                                 (Handle(GEOM_Object)              theMainShape,
1469                                  Handle(TColStd_HArray1OfInteger) theIndices)
1470 {
1471   SetErrorCode(KO);
1472
1473   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
1474
1475   if (!theIndices->Length()) {
1476     SetErrorCode(NOT_FOUND_ANY);
1477     return aSeq;
1478   }
1479
1480   if (theMainShape.IsNull()) return NULL;
1481   TopoDS_Shape aShape = theMainShape->GetValue();
1482   if (aShape.IsNull()) return NULL;
1483
1484   Handle(GEOM_Function) aMainShape = theMainShape->GetLastFunction();
1485
1486   TopTools_IndexedMapOfShape anIndices;
1487   TopExp::MapShapes(aShape, anIndices);
1488
1489   Handle(TColStd_HArray1OfInteger) anArray;
1490   Handle(GEOM_Object) anObj;
1491
1492   TCollection_AsciiString anAsciiList, anEntry;
1493   Standard_Integer i, low = theIndices->Lower(), up = theIndices->Upper();
1494   for (i = low; i <= up; i++) {
1495     int id = theIndices->Value(i);
1496     if (1 <= id && id <= anIndices.Extent()) {
1497       TopoDS_Shape aValue = anIndices.FindKey(id);
1498       anArray = new TColStd_HArray1OfInteger(1,1);
1499       anArray->SetValue(1, id);
1500
1501       anObj = GetEngine()->AddObject(GetDocID(), GEOM_SUBSHAPE);
1502       if (!anObj.IsNull()) {
1503         Handle(GEOM_Function) aFunction = anObj->AddFunction(GEOM_Object::GetSubShapeID(), 1);
1504         if (aFunction.IsNull()) return aSeq;
1505
1506         GEOM_ISubShape aSSI (aFunction);
1507         aSSI.SetMainShape(aMainShape);
1508         aSSI.SetIndices(anArray);
1509
1510         // Set function value directly, as we know it.
1511         // Usage of Solver here would lead to significant loss of time,
1512         // because GEOM_SubShapeDriver will build TopTools_IndexedMapOfShape
1513         // on the main shape for each being calculated sub-shape separately.
1514         aFunction->SetValue(aValue);
1515
1516         // Put this sub-shape in the list of sub-shapes of theMainShape
1517         aMainShape->AddSubShapeReference(aFunction);
1518
1519         aSeq->Append(anObj);
1520
1521         // for python command
1522         TDF_Tool::Entry(anObj->GetEntry(), anEntry);
1523         anAsciiList += anEntry;
1524         anAsciiList += ",";
1525       }
1526     }
1527   }
1528
1529   //Make a Python command
1530   anAsciiList.Trunc(anAsciiList.Length() - 1);
1531
1532   GEOM::TPythonDump pd (aMainShape, /*append=*/true);
1533   pd << "[" << anAsciiList.ToCString() << "] = geompy.SubShapes("
1534      << theMainShape << ", [" ;
1535   for (i = low; i <= up - 1; i++) {
1536     pd << theIndices->Value(i) << ", ";
1537   }
1538   pd << theIndices->Value(up) << "])";
1539
1540   SetErrorCode(OK);
1541
1542   return aSeq;
1543 }
1544
1545 //=============================================================================
1546 /*!
1547  *  GetSubShapeIndex
1548  */
1549 //=============================================================================
1550 Standard_Integer GEOMImpl_IShapesOperations::GetSubShapeIndex (Handle(GEOM_Object) theMainShape,
1551                                                                Handle(GEOM_Object) theSubShape)
1552 {
1553   SetErrorCode(KO);
1554
1555   TopoDS_Shape aMainShape = theMainShape->GetValue();
1556   TopoDS_Shape aSubShape = theSubShape->GetValue();
1557
1558   if (aMainShape.IsNull() || aSubShape.IsNull()) return -1;
1559
1560   TopTools_IndexedMapOfShape anIndices;
1561   TopExp::MapShapes(aMainShape, anIndices);
1562   if (anIndices.Contains(aSubShape)) {
1563     SetErrorCode(OK);
1564     return anIndices.FindIndex(aSubShape);
1565   }
1566
1567   return -1;
1568 }
1569
1570 //=============================================================================
1571 /*!
1572  *  GetTopologyIndex
1573  */
1574 //=============================================================================
1575 Standard_Integer GEOMImpl_IShapesOperations::GetTopologyIndex (Handle(GEOM_Object) theMainShape,
1576                                                                Handle(GEOM_Object) theSubShape)
1577 {
1578   SetErrorCode(OK);
1579
1580   TopoDS_Shape aMainShape = theMainShape->GetValue();
1581   TopoDS_Shape aSubShape = theSubShape->GetValue();
1582
1583   if (aMainShape.IsNull() || aSubShape.IsNull()) {
1584     SetErrorCode("Null argument shape given");
1585     return -1;
1586   }
1587
1588   int index = 1;
1589   if (aSubShape.ShapeType() == TopAbs_COMPOUND) {
1590     TopoDS_Iterator it;
1591     TopTools_ListOfShape CL;
1592     CL.Append(aMainShape);
1593     TopTools_ListIteratorOfListOfShape itC;
1594     for (itC.Initialize(CL); itC.More(); itC.Next()) {
1595       for (it.Initialize(itC.Value()); it.More(); it.Next()) {
1596         if (it.Value().ShapeType() == TopAbs_COMPOUND) {
1597           if (it.Value().IsSame(aSubShape))
1598             return index;
1599           else
1600             index++;
1601           CL.Append(it.Value());
1602         }
1603       }
1604     }
1605   } else {
1606     TopExp_Explorer anExp (aMainShape, aSubShape.ShapeType());
1607     TopTools_MapOfShape M;
1608     for (; anExp.More(); anExp.Next()) {
1609       if (M.Add(anExp.Current())) {
1610         if (anExp.Current().IsSame(aSubShape))
1611           return index;
1612         index++;
1613       }
1614     }
1615   }
1616
1617   SetErrorCode("The sub-shape does not belong to the main shape");
1618   return -1;
1619 }
1620
1621 //=============================================================================
1622 /*!
1623  *  GetShapeTypeString
1624  */
1625 //=============================================================================
1626 TCollection_AsciiString GEOMImpl_IShapesOperations::GetShapeTypeString (Handle(GEOM_Object) theShape)
1627 {
1628   SetErrorCode(KO);
1629
1630   TCollection_AsciiString aTypeName ("Null Shape");
1631
1632   TopoDS_Shape aShape = theShape->GetValue();
1633   if (aShape.IsNull())
1634     return aTypeName;
1635
1636   switch (aShape.ShapeType() )
1637   {
1638   case TopAbs_COMPOUND:
1639     aTypeName = "Compound";
1640     break;
1641   case  TopAbs_COMPSOLID:
1642     aTypeName = "Compound Solid";
1643     break;
1644   case TopAbs_SOLID:
1645     aTypeName = "Solid";
1646     break;
1647   case TopAbs_SHELL:
1648     aTypeName = "Shell";
1649     break;
1650   case TopAbs_FACE:
1651     {
1652       BRepAdaptor_Surface surf (TopoDS::Face(aShape));
1653       if (surf.GetType() == GeomAbs_Plane)
1654         aTypeName = "Plane";
1655       else if (surf.GetType() == GeomAbs_Cylinder)
1656         aTypeName = "Cylindrical Face";
1657       else if (surf.GetType() == GeomAbs_Sphere)
1658         aTypeName = "Spherical Face";
1659       else if (surf.GetType() == GeomAbs_Torus)
1660         aTypeName = "Toroidal Face";
1661       else if (surf.GetType() == GeomAbs_Cone)
1662         aTypeName = "Conical Face";
1663       else
1664         aTypeName = "GEOM::FACE";
1665     }
1666     break;
1667   case TopAbs_WIRE:
1668     aTypeName = "Wire";
1669     break;
1670   case TopAbs_EDGE:
1671     {
1672       BRepAdaptor_Curve curv (TopoDS::Edge(aShape));
1673       if (curv.GetType() == GeomAbs_Line) {
1674         if ((Abs(curv.FirstParameter()) >= 1E6) ||
1675             (Abs(curv.LastParameter()) >= 1E6))
1676           aTypeName = "Line";
1677         else
1678           aTypeName = "Edge";
1679       } else if (curv.GetType() == GeomAbs_Circle) {
1680         if (curv.IsClosed())
1681           aTypeName = "Circle";
1682         else
1683           aTypeName = "Arc";
1684       } else {
1685         aTypeName = "Edge";
1686       }
1687     }
1688     break;
1689   case TopAbs_VERTEX:
1690     aTypeName = "Vertex";
1691     break;
1692   case TopAbs_SHAPE:
1693     aTypeName = "Shape";
1694     break;
1695   default:
1696     aTypeName = "Shape of unknown type";
1697   }
1698
1699   return aTypeName;
1700 }
1701
1702 //=============================================================================
1703 /*!
1704  *  NumberOfSubShapes
1705  */
1706 //=============================================================================
1707 Standard_Integer GEOMImpl_IShapesOperations::NumberOfSubShapes
1708                                           (Handle(GEOM_Object)    theShape,
1709                                            const Standard_Integer theShapeType)
1710 {
1711   SetErrorCode(KO);
1712   Standard_Integer nbShapes = 0;
1713
1714   if (theShape.IsNull()) return -1;
1715   TopoDS_Shape aShape = theShape->GetValue();
1716   if (aShape.IsNull()) return -1;
1717
1718   /*
1719   TopTools_MapOfShape mapShape;
1720
1721   if (aShape.ShapeType() == TopAbs_COMPOUND &&
1722       (TopAbs_ShapeEnum(theShapeType) == TopAbs_SHAPE ||
1723        TopAbs_ShapeEnum(theShapeType) == TopAbs_COMPSOLID ||
1724        TopAbs_ShapeEnum(theShapeType) == TopAbs_COMPOUND)) {
1725     TopoDS_Iterator It (aShape, Standard_True, Standard_True);
1726     for (; It.More(); It.Next()) {
1727       if (mapShape.Add(It.Value())) {
1728         if (TopAbs_ShapeEnum(theShapeType) == TopAbs_SHAPE ||
1729             TopAbs_ShapeEnum(theShapeType) == It.Value().ShapeType()) {
1730           nbShapes++;
1731         }
1732       }
1733     }
1734   } else {
1735     TopExp_Explorer exp (aShape, TopAbs_ShapeEnum(theShapeType));
1736     for (; exp.More(); exp.Next())
1737       if (mapShape.Add(exp.Current()))
1738         nbShapes++;
1739   }
1740   */
1741
1742   try {
1743 #if OCC_VERSION_LARGE > 0x06010000
1744     OCC_CATCH_SIGNALS;
1745 #endif
1746     int iType, nbTypes [TopAbs_SHAPE];
1747     for (iType = 0; iType < TopAbs_SHAPE; ++iType)
1748       nbTypes[iType] = 0;
1749     nbTypes[aShape.ShapeType()]++;
1750
1751     TopTools_MapOfShape aMapOfShape;
1752     aMapOfShape.Add(aShape);
1753     TopTools_ListOfShape aListOfShape;
1754     aListOfShape.Append(aShape);
1755
1756     TopTools_ListIteratorOfListOfShape itL (aListOfShape);
1757     for (; itL.More(); itL.Next()) {
1758       TopoDS_Iterator it (itL.Value());
1759       for (; it.More(); it.Next()) {
1760         TopoDS_Shape s = it.Value();
1761         if (aMapOfShape.Add(s)) {
1762           aListOfShape.Append(s);
1763           nbTypes[s.ShapeType()]++;
1764         }
1765       }
1766     }
1767
1768     if (TopAbs_ShapeEnum(theShapeType) == TopAbs_SHAPE)
1769       nbShapes = aMapOfShape.Extent();
1770     else
1771       nbShapes = nbTypes[theShapeType];
1772   }
1773   catch (Standard_Failure) {
1774     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1775     SetErrorCode(aFail->GetMessageString());
1776     return -1;
1777   }
1778
1779   SetErrorCode(OK);
1780   return nbShapes;
1781 }
1782
1783 //=============================================================================
1784 /*!
1785  *  ReverseShape
1786  */
1787 //=============================================================================
1788 Handle(GEOM_Object) GEOMImpl_IShapesOperations::ReverseShape(Handle(GEOM_Object) theShape)
1789 {
1790   SetErrorCode(KO);
1791
1792   if (theShape.IsNull()) return NULL;
1793
1794   /*
1795   //Add a new reversed object
1796   Handle(GEOM_Object) aReversed = GetEngine()->AddObject(GetDocID(), theShape->GetType());
1797
1798   //Add a new Revese function
1799   Handle(GEOM_Function) aFunction;
1800   aFunction = aReversed->AddFunction(GEOMImpl_ShapeDriver::GetID(), REVERSE_ORIENTATION);
1801   if (aFunction.IsNull()) return NULL;
1802
1803   //Check if the function is set correctly
1804   if (aFunction->GetDriverGUID() != GEOMImpl_ShapeDriver::GetID()) return NULL;
1805
1806   GEOMImpl_IShapes aSI (aFunction);
1807
1808   Handle(GEOM_Function) aRefShape = theShape->GetLastFunction();
1809   if (aRefShape.IsNull()) return NULL;
1810
1811   aSI.SetBase(aRefShape);
1812
1813   //Compute the sub-shape value
1814   try {
1815 #if OCC_VERSION_LARGE > 0x06010000
1816     OCC_CATCH_SIGNALS;
1817 #endif
1818     if (!GetSolver()->ComputeFunction(aFunction)) {
1819       SetErrorCode("Shape driver failed to reverse shape");
1820       return NULL;
1821     }
1822   }
1823   catch (Standard_Failure) {
1824     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1825     SetErrorCode(aFail->GetMessageString());
1826     return NULL;
1827   }
1828
1829   //Make a Python command
1830   GEOM::TPythonDump(aFunction) << aReversed
1831     << " = geompy.ChangeOrientation(" << theShape << ")";
1832
1833   SetErrorCode(OK);
1834   */
1835
1836   Handle(GEOM_Object) aReversed;
1837
1838   GEOM_Engine* anEngine = GetEngine();
1839   //GEOMImpl_Gen* aGen = dynamic_cast<GEOMImpl_Gen*>(anEngine);
1840   GEOMImpl_Gen* aGen = (GEOMImpl_Gen*)anEngine;
1841
1842   if (aGen) {
1843     GEOMImpl_IHealingOperations* anIHealingOperations =
1844       aGen->GetIHealingOperations(GetDocID());
1845     aReversed = anIHealingOperations->ChangeOrientationCopy(theShape);
1846     SetErrorCode(anIHealingOperations->GetErrorCode());
1847   }
1848
1849   return aReversed;
1850 }
1851
1852 //=============================================================================
1853 /*!
1854  *  GetFreeFacesIDs
1855  */
1856 //=============================================================================
1857 Handle(TColStd_HSequenceOfInteger) GEOMImpl_IShapesOperations::GetFreeFacesIDs
1858                                                  (Handle(GEOM_Object) theShape)
1859 {
1860   SetErrorCode(KO);
1861
1862   if (theShape.IsNull()) return NULL;
1863   TopoDS_Shape aShape = theShape->GetValue();
1864   if (aShape.IsNull()) return NULL;
1865
1866   Handle(TColStd_HSequenceOfInteger) aSeq = new TColStd_HSequenceOfInteger;
1867
1868   TopTools_IndexedDataMapOfShapeListOfShape mapFaceBlocks;
1869   GEOMImpl_Block6Explorer::MapShapesAndAncestors
1870     (aShape, TopAbs_FACE, TopAbs_SOLID, mapFaceBlocks);
1871
1872   Standard_Integer ind = 1, nbFaces = mapFaceBlocks.Extent();
1873
1874   if (nbFaces == 0) {
1875     SetErrorCode("The given shape has no faces");
1876     return aSeq;
1877   }
1878
1879   TopTools_IndexedMapOfShape anIndices;
1880   TopExp::MapShapes(aShape, anIndices);
1881
1882   Standard_Integer id;
1883   for (; ind <= nbFaces; ind++) {
1884     if (mapFaceBlocks.FindFromIndex(ind).Extent() != 2) {
1885       id = anIndices.FindIndex(mapFaceBlocks.FindKey(ind));
1886       aSeq->Append(id);
1887     }
1888   }
1889
1890   //The explode doesn't change object so no new function is required.
1891   Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
1892
1893   //Make a Python command
1894   GEOM::TPythonDump(aFunction, /*append=*/true)
1895     << "listFreeFacesIDs = geompy.GetFreeFacesIDs(" << theShape << ")";
1896
1897   SetErrorCode(OK);
1898   return aSeq;
1899 }
1900
1901 //=======================================================================
1902 //function : GetSharedShapes
1903 //purpose  :
1904 //=======================================================================
1905 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetSharedShapes
1906                                                 (Handle(GEOM_Object)    theShape1,
1907                                                  Handle(GEOM_Object)    theShape2,
1908                                                  const Standard_Integer theShapeType)
1909 {
1910   SetErrorCode(KO);
1911
1912   if (theShape1.IsNull() || theShape2.IsNull()) return NULL;
1913
1914   TopoDS_Shape aShape1 = theShape1->GetValue();
1915   TopoDS_Shape aShape2 = theShape2->GetValue();
1916
1917   if (aShape1.IsNull() || aShape2.IsNull()) return NULL;
1918
1919   TopTools_IndexedMapOfShape anIndices;
1920   TopExp::MapShapes(aShape1, anIndices);
1921   Handle(TColStd_HArray1OfInteger) anArray;
1922
1923   TopTools_IndexedMapOfShape mapShape1;
1924   TopExp::MapShapes(aShape1, TopAbs_ShapeEnum(theShapeType), mapShape1);
1925
1926   Handle(GEOM_Object) anObj;
1927   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
1928   TCollection_AsciiString anAsciiList, anEntry;
1929
1930   TopTools_MapOfShape mapShape2;
1931   TopExp_Explorer exp (aShape2, TopAbs_ShapeEnum(theShapeType));
1932   for (; exp.More(); exp.Next()) {
1933     TopoDS_Shape aSS = exp.Current();
1934     if (mapShape2.Add(aSS) && mapShape1.Contains(aSS)) {
1935       anArray = new TColStd_HArray1OfInteger(1,1);
1936       anArray->SetValue(1, anIndices.FindIndex(aSS));
1937       anObj = GetEngine()->AddSubShape(theShape1, anArray);
1938       aSeq->Append(anObj);
1939
1940       // for python command
1941       TDF_Tool::Entry(anObj->GetEntry(), anEntry);
1942       anAsciiList += anEntry;
1943       anAsciiList += ",";
1944     }
1945   }
1946
1947   if (aSeq->IsEmpty()) {
1948     SetErrorCode("The given shapes have no shared sub-shapes of the requested type");
1949     return aSeq;
1950   }
1951
1952   //Make a Python command
1953   anAsciiList.Trunc(anAsciiList.Length() - 1);
1954
1955   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
1956
1957   GEOM::TPythonDump(aFunction) << "[" << anAsciiList.ToCString()
1958     << "] = geompy.GetSharedShapes(" << theShape1 << ", "
1959       << theShape2 << ", " << TopAbs_ShapeEnum(theShapeType) << ")";
1960
1961   SetErrorCode(OK);
1962   return aSeq;
1963 }
1964
1965 //=======================================================================
1966 //function : GetSharedShapes
1967 //purpose  :
1968 //=======================================================================
1969 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetSharedShapes
1970                                      (std::list<Handle(GEOM_Object)> theShapes,
1971                                       const Standard_Integer         theShapeType)
1972 {
1973   SetErrorCode(KO);
1974
1975   int aLen = theShapes.size();
1976   if (aLen < 1) return NULL;
1977
1978   int ind = 1;
1979   std::list<Handle(GEOM_Object)>::iterator it = theShapes.begin();
1980
1981   Handle(GEOM_Object) aMainObj = (*it++);
1982   Handle(GEOM_Function) aMainShape = aMainObj->GetLastFunction();
1983   if (aMainShape.IsNull()) {
1984     SetErrorCode("NULL shape for GetSharedShapes");
1985     return NULL;
1986   }
1987
1988   TopoDS_Shape aShape1 = aMainShape->GetValue();
1989   if (aShape1.IsNull()) return NULL;
1990
1991   TopTools_IndexedMapOfShape anIndices;
1992   TopExp::MapShapes(aShape1, anIndices);
1993
1994   TopTools_IndexedMapOfShape mapSelected;
1995   TopExp::MapShapes(aShape1, TopAbs_ShapeEnum(theShapeType), mapSelected);
1996
1997   // Find shared shapes
1998   BRep_Builder B;
1999   TopoDS_Compound aCurrSelection;
2000
2001   for (; it != theShapes.end(); it++, ind++) {
2002     Handle(GEOM_Function) aRefShape = (*it)->GetLastFunction();
2003     if (aRefShape.IsNull()) {
2004       SetErrorCode("NULL shape for GetSharedShapes");
2005       return NULL;
2006     }
2007
2008     TopoDS_Compound aCompound;
2009     B.MakeCompound(aCompound);
2010
2011     TopoDS_Shape aShape2 = aRefShape->GetValue();
2012     if (aShape2.IsNull()) return NULL;
2013
2014     TopTools_MapOfShape mapShape2;
2015     TopExp_Explorer exp (aShape2, TopAbs_ShapeEnum(theShapeType));
2016     for (; exp.More(); exp.Next()) {
2017       TopoDS_Shape aSS = exp.Current();
2018       if (mapShape2.Add(aSS) && mapSelected.Contains(aSS)) {
2019         B.Add(aCompound, aSS);
2020       }
2021     }
2022
2023     mapSelected.Clear();
2024     TopExp::MapShapes(aCompound, TopAbs_ShapeEnum(theShapeType), mapSelected);
2025     aCurrSelection = aCompound;
2026   }
2027
2028   // Create GEOM_Object for each found shared shape (collected in aCurrSelection)
2029   Handle(GEOM_Object) anObj;
2030   Handle(TColStd_HArray1OfInteger) anArray;
2031   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
2032   TCollection_AsciiString anAsciiList, anEntry;
2033
2034   TopoDS_Iterator itSel (aCurrSelection, Standard_True, Standard_True);
2035   for (; itSel.More(); itSel.Next()) {
2036     anArray = new TColStd_HArray1OfInteger(1,1);
2037     anArray->SetValue(1, anIndices.FindIndex(itSel.Value()));
2038     anObj = GetEngine()->AddSubShape(aMainObj, anArray);
2039     aSeq->Append(anObj);
2040
2041     // for python command
2042     TDF_Tool::Entry(anObj->GetEntry(), anEntry);
2043     anAsciiList += anEntry;
2044     anAsciiList += ",";
2045   }
2046
2047   if (aSeq->IsEmpty()) {
2048     SetErrorCode("The given shapes have no shared sub-shapes of the requested type");
2049     return aSeq;
2050   }
2051
2052   // Make a Python command
2053   anAsciiList.Trunc(anAsciiList.Length() - 1);
2054
2055   // IPAL22904: TC6.5.0: order of python commands is wrong after dump study
2056   Handle(TColStd_HSequenceOfTransient) anObjects = new TColStd_HSequenceOfTransient;
2057   for( it = theShapes.begin(); it != theShapes.end(); it++ )
2058   {
2059     Handle(GEOM_Object) anObj = *it;
2060     if( !anObj.IsNull() )
2061       anObjects->Append( anObj );
2062   }
2063
2064   // Get the function of the latest published object
2065   Handle(GEOM_Function) aFunction = GEOM::GetCreatedLast( anObjects )->GetLastFunction();
2066   if( aFunction.IsNull() ) // just in case
2067     aFunction = aMainShape;
2068
2069   GEOM::TPythonDump pd (aFunction, /*append=*/true);
2070   pd << "[" << anAsciiList.ToCString()
2071      << "] = geompy.GetSharedShapesMulti([";
2072
2073   it = theShapes.begin();
2074   pd << (*it++);
2075   while (it != theShapes.end()) {
2076     pd << ", " << (*it++);
2077   }
2078
2079   pd << "], " << TopAbs_ShapeEnum(theShapeType) << ")";
2080
2081   SetErrorCode(OK);
2082   return aSeq;
2083 }
2084
2085 //=============================================================================
2086 /*!
2087  *
2088  */
2089 //=============================================================================
2090 static GEOM::TPythonDump& operator<< (GEOM::TPythonDump&   theDump,
2091                                       const GEOMAlgo_State theState)
2092 {
2093   switch (theState) {
2094   case GEOMAlgo_ST_IN:
2095     theDump << "geompy.GEOM.ST_IN";
2096     break;
2097   case GEOMAlgo_ST_OUT:
2098     theDump << "geompy.GEOM.ST_OUT";
2099     break;
2100   case GEOMAlgo_ST_ON:
2101     theDump << "geompy.GEOM.ST_ON";
2102     break;
2103   case GEOMAlgo_ST_ONIN:
2104     theDump << "geompy.GEOM.ST_ONIN";
2105     break;
2106   case GEOMAlgo_ST_ONOUT:
2107     theDump << "geompy.GEOM.ST_ONOUT";
2108     break;
2109   default:
2110     theDump << "geompy.GEOM.ST_UNKNOWN";
2111     break;
2112   }
2113   return theDump;
2114 }
2115
2116 //=======================================================================
2117 //function : checkTypeShapesOn
2118 /*!
2119  * \brief Checks if theShapeType parameter of GetShapesOnXXX() is OK
2120  * \param theShapeType - the shape type to check
2121  * \retval bool  - result of the check
2122  */
2123 //=======================================================================
2124 bool GEOMImpl_IShapesOperations::checkTypeShapesOn(const Standard_Integer theShapeType)
2125 {
2126   if (theShapeType != TopAbs_VERTEX &&
2127       theShapeType != TopAbs_EDGE &&
2128       theShapeType != TopAbs_FACE &&
2129       theShapeType != TopAbs_SOLID) {
2130     SetErrorCode("Only solids, vertices, edges or faces can be found by this method");
2131     return false;
2132   }
2133   return true;
2134 }
2135
2136 //=======================================================================
2137 //function : makePlane
2138   /*!
2139    * \brief Creates Geom_Plane
2140     * \param theAx1 - shape object defining plane parameters
2141     * \retval Handle(Geom_Surface) - resulting surface
2142    */
2143 //=======================================================================
2144 Handle(Geom_Surface) GEOMImpl_IShapesOperations::makePlane(const TopoDS_Shape& anAx1)
2145 {
2146   if (anAx1.ShapeType() != TopAbs_EDGE) return NULL;
2147   TopoDS_Edge anEdge = TopoDS::Edge(anAx1);
2148   TopoDS_Vertex V1, V2;
2149   TopExp::Vertices(anEdge, V1, V2, Standard_True);
2150   if (V1.IsNull() || V2.IsNull()) {
2151     SetErrorCode("Bad edge given for the plane normal vector");
2152     return NULL;
2153   }
2154   gp_Pnt aLoc = BRep_Tool::Pnt(V1);
2155   gp_Vec aVec (aLoc, BRep_Tool::Pnt(V2));
2156   if (aVec.Magnitude() < Precision::Confusion()) {
2157     SetErrorCode("Vector with null magnitude given");
2158     return NULL;
2159   }
2160   return new Geom_Plane(aLoc, aVec);
2161 }
2162
2163 //=======================================================================
2164 //function : makeCylinder
2165   /*!
2166    * \brief Creates Geom_CylindricalSurface
2167     * \param theAx1 - edge defining cylinder axis
2168     * \param theRadius - cylinder radius
2169     * \retval Handle(Geom_Surface) - resulting surface
2170    */
2171 //=======================================================================
2172 Handle(Geom_Surface) GEOMImpl_IShapesOperations::makeCylinder(const TopoDS_Shape& anAxis,
2173                                                               const Standard_Real theRadius)
2174 {
2175   //Axis of the cylinder
2176   if (anAxis.ShapeType() != TopAbs_EDGE) {
2177     SetErrorCode("Not an edge given for the axis");
2178     return NULL;
2179   }
2180   TopoDS_Edge anEdge = TopoDS::Edge(anAxis);
2181   TopoDS_Vertex V1, V2;
2182   TopExp::Vertices(anEdge, V1, V2, Standard_True);
2183   if (V1.IsNull() || V2.IsNull()) {
2184     SetErrorCode("Bad edge given for the axis");
2185     return NULL;
2186   }
2187   gp_Pnt aLoc = BRep_Tool::Pnt(V1);
2188   gp_Vec aVec (aLoc, BRep_Tool::Pnt(V2));
2189   if (aVec.Magnitude() < Precision::Confusion()) {
2190     SetErrorCode("Vector with null magnitude given");
2191     return NULL;
2192   }
2193
2194   gp_Ax3 anAx3 (aLoc, aVec);
2195   return new Geom_CylindricalSurface(anAx3, theRadius);
2196 }
2197
2198 //=======================================================================
2199 //function : getShapesOnBoxIDs
2200   /*!
2201    * \brief Find IDs of sub-shapes complying with given status about surface
2202     * \param theBox - the box to check state of sub-shapes against
2203     * \param theShape - the shape to explore
2204     * \param theShapeType - type of sub-shape of theShape
2205     * \param theState - required state
2206     * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found sub-shapes
2207    */
2208 //=======================================================================
2209 Handle(TColStd_HSequenceOfInteger)
2210   GEOMImpl_IShapesOperations::getShapesOnBoxIDs(const Handle(GEOM_Object)& theBox,
2211                                                 const Handle(GEOM_Object)& theShape,
2212                                                 const Standard_Integer theShapeType,
2213                                                 GEOMAlgo_State theState)
2214 {
2215   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs;
2216
2217   TopoDS_Shape aBox = theBox->GetValue();
2218   TopoDS_Shape aShape = theShape->GetValue();
2219
2220   // Check presence of triangulation, build if need
2221   if (!CheckTriangulation(aShape)) {
2222     SetErrorCode("Cannot build triangulation on the shape");
2223     return aSeqOfIDs;
2224   }
2225
2226   // Call algo
2227   GEOMAlgo_FinderShapeOn2 aFinder;
2228   Standard_Real aTol = 0.0001; // default value
2229
2230   Handle(GEOMAlgo_ClsfBox) aClsfBox = new GEOMAlgo_ClsfBox;
2231   aClsfBox->SetBox(aBox);
2232
2233   aFinder.SetShape(aShape);
2234   aFinder.SetTolerance(aTol);
2235   aFinder.SetClsf(aClsfBox);
2236   aFinder.SetShapeType( (TopAbs_ShapeEnum)theShapeType );
2237   aFinder.SetState(theState);
2238   aFinder.Perform();
2239
2240   // Interprete results
2241   Standard_Integer iErr = aFinder.ErrorStatus();
2242   // the detailed description of error codes is in GEOMAlgo_FinderShapeOn1.cxx
2243   if (iErr) {
2244     MESSAGE(" iErr : " << iErr);
2245     TCollection_AsciiString aMsg (" iErr : ");
2246     aMsg += TCollection_AsciiString(iErr);
2247     SetErrorCode(aMsg);
2248     return aSeqOfIDs;
2249   }
2250   Standard_Integer iWrn = aFinder.WarningStatus();
2251   // the detailed description of warning codes is in GEOMAlgo_FinderShapeOn1.cxx
2252   if (iWrn) {
2253     MESSAGE(" *** iWrn : " << iWrn);
2254   }
2255
2256   const TopTools_ListOfShape& listSS = aFinder.Shapes(); // the result
2257
2258   if (listSS.Extent() < 1) {
2259     //SetErrorCode("Not a single sub-shape of the requested type found on the given surface");
2260     SetErrorCode(NOT_FOUND_ANY); // NPAL18017
2261     return aSeqOfIDs;
2262   }
2263
2264   // Fill sequence of object IDs
2265   aSeqOfIDs = new TColStd_HSequenceOfInteger;
2266
2267   TopTools_IndexedMapOfShape anIndices;
2268   TopExp::MapShapes(aShape, anIndices);
2269
2270   TopTools_ListIteratorOfListOfShape itSub (listSS);
2271   for (int index = 1; itSub.More(); itSub.Next(), ++index) {
2272     int id = anIndices.FindIndex(itSub.Value());
2273     aSeqOfIDs->Append(id);
2274   }
2275
2276   return aSeqOfIDs;
2277 }
2278
2279 //=======================================================================
2280 //function : GetShapesOnBoxIDs
2281 /*!
2282    * \brief Find sub-shapes complying with given status about surface
2283     * \param theBox - the box to check state of sub-shapes against
2284     * \param theShape - the shape to explore
2285     * \param theShapeType - type of sub-shape of theShape
2286     * \param theState - required state
2287     * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found sub-shapes
2288  */
2289 //=======================================================================
2290 Handle(TColStd_HSequenceOfInteger)
2291     GEOMImpl_IShapesOperations::GetShapesOnBoxIDs(const Handle(GEOM_Object)& theBox,
2292                                                   const Handle(GEOM_Object)& theShape,
2293                                                   const Standard_Integer theShapeType,
2294                                                   GEOMAlgo_State theState)
2295 {
2296   // Find sub-shapes ids
2297   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
2298     getShapesOnBoxIDs (theBox, theShape, theShapeType, theState);
2299   if ( aSeqOfIDs.IsNull()  || aSeqOfIDs->Length() == 0 )
2300     return NULL;
2301
2302   // The GetShapesOnBox() doesn't change object so no new function is required.
2303   Handle(GEOM_Function) aFunction = GEOM::GetCreatedLast(theShape,theBox)->GetLastFunction();
2304
2305   // Make a Python command
2306   GEOM::TPythonDump(aFunction)
2307     << "listShapesOnBoxIDs = geompy.GetShapesOnBoxIDs("
2308     << theBox << ", "
2309     << theShape << ", "
2310     << TopAbs_ShapeEnum(theShapeType) << ", "
2311     << theState << ")";
2312
2313   SetErrorCode(OK);
2314   return aSeqOfIDs;
2315 }
2316
2317 //=======================================================================
2318 //function : GetShapesOnBox
2319 /*!
2320    * \brief Find sub-shapes complying with given status about surface
2321     * \param theBox - the box to check state of sub-shapes against
2322     * \param theShape - the shape to explore
2323     * \param theShapeType - type of sub-shape of theShape
2324     * \param theState - required state
2325     * \retval Handle(TColStd_HSequenceOfTransient) - found sub-shapes
2326  */
2327 //=======================================================================
2328 Handle(TColStd_HSequenceOfTransient)
2329     GEOMImpl_IShapesOperations::GetShapesOnBox(const Handle(GEOM_Object)& theBox,
2330                                                const Handle(GEOM_Object)&  theShape,
2331                                                const Standard_Integer theShapeType,
2332                                                GEOMAlgo_State theState)
2333 {
2334   // Find sub-shapes ids
2335   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
2336     getShapesOnBoxIDs (theBox, theShape, theShapeType, theState);
2337   if ( aSeqOfIDs.IsNull()  || aSeqOfIDs->Length() == 0 )
2338     return NULL;
2339
2340   // Find objects by indices
2341   TCollection_AsciiString anAsciiList;
2342   Handle(TColStd_HSequenceOfTransient) aSeq;
2343   aSeq = getObjectsShapesOn( theShape, aSeqOfIDs, anAsciiList );
2344   if ( aSeq.IsNull() || aSeq->IsEmpty() )
2345     return NULL;
2346
2347   // Make a Python command
2348
2349   Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast( aSeq->Value( 1 ));
2350   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
2351
2352   GEOM::TPythonDump(aFunction)
2353     << "[" << anAsciiList.ToCString() << "] = geompy.GetShapesOnBox("
2354     << theBox << ", "
2355     << theShape << ", "
2356     << TopAbs_ShapeEnum(theShapeType) << ", "
2357     << theState << ")";
2358
2359   SetErrorCode(OK);
2360   return aSeq;
2361 }
2362
2363 //=======================================================================
2364 //function : getShapesOnShapeIDs
2365 /*!
2366  * \brief Find IDs of sub-shapes complying with given status about surface
2367  * \param theCheckShape - the shape to check state of sub-shapes against
2368  * \param theShape - the shape to explore
2369  * \param theShapeType - type of sub-shape of theShape
2370  * \param theState - required state
2371  * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found sub-shapes
2372  */
2373 //=======================================================================
2374 Handle(TColStd_HSequenceOfInteger)
2375   GEOMImpl_IShapesOperations::getShapesOnShapeIDs
2376                                  (const Handle(GEOM_Object)& theCheckShape,
2377                                   const Handle(GEOM_Object)& theShape,
2378                                   const Standard_Integer theShapeType,
2379                                   GEOMAlgo_State theState)
2380 {
2381   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs;
2382
2383   TopoDS_Shape aCheckShape = theCheckShape->GetValue();
2384   TopoDS_Shape aShape = theShape->GetValue();
2385   TopTools_ListOfShape res;
2386
2387   // Check presence of triangulation, build if need
2388   if (!CheckTriangulation(aShape)) {
2389     SetErrorCode("Cannot build triangulation on the shape");
2390     return aSeqOfIDs;
2391   }
2392
2393   // Call algo
2394   GEOMAlgo_FinderShapeOn2 aFinder;
2395   Standard_Real aTol = 0.0001; // default value
2396
2397   Handle(GEOMAlgo_ClsfSolid) aClsfSolid = new GEOMAlgo_ClsfSolid;
2398   aClsfSolid->SetShape(aCheckShape);
2399
2400   aFinder.SetShape(aShape);
2401   aFinder.SetTolerance(aTol);
2402   aFinder.SetClsf(aClsfSolid);
2403   aFinder.SetShapeType( (TopAbs_ShapeEnum)theShapeType );
2404   aFinder.SetState(theState);
2405   aFinder.Perform();
2406
2407   // Interprete results
2408   Standard_Integer iErr = aFinder.ErrorStatus();
2409   // the detailed description of error codes is in GEOMAlgo_FinderShapeOn1.cxx
2410   if (iErr) {
2411     if (iErr == 41) {
2412       SetErrorCode("theCheckShape must be a solid");
2413     }
2414     else {
2415       MESSAGE(" iErr : " << iErr);
2416       TCollection_AsciiString aMsg (" iErr : ");
2417       aMsg += TCollection_AsciiString(iErr);
2418       SetErrorCode(aMsg);
2419     }
2420     return aSeqOfIDs;
2421   }
2422   Standard_Integer iWrn = aFinder.WarningStatus();
2423   // the detailed description of warning codes is in GEOMAlgo_FinderShapeOn1.cxx
2424   if (iWrn) {
2425     MESSAGE(" *** iWrn : " << iWrn);
2426   }
2427
2428   const TopTools_ListOfShape& listSS = aFinder.Shapes(); // the result
2429
2430   if (listSS.Extent() < 1) {
2431     //SetErrorCode("Not a single sub-shape of the requested type found on the given surface");
2432     SetErrorCode(NOT_FOUND_ANY); // NPAL18017
2433   }
2434
2435   // Fill sequence of object IDs
2436   aSeqOfIDs = new TColStd_HSequenceOfInteger;
2437
2438   TopTools_IndexedMapOfShape anIndices;
2439   TopExp::MapShapes(aShape, anIndices);
2440
2441   TopTools_ListIteratorOfListOfShape itSub (listSS);
2442   for (int index = 1; itSub.More(); itSub.Next(), ++index) {
2443     int id = anIndices.FindIndex(itSub.Value());
2444     aSeqOfIDs->Append(id);
2445   }
2446
2447   return aSeqOfIDs;
2448 }
2449
2450 //=======================================================================
2451 //function : GetShapesOnShapeIDs
2452 /*!
2453  * \brief Find sub-shapes complying with given status about surface
2454  * \param theCheckShape - the shape to check state of sub-shapes against
2455  * \param theShape - the shape to explore
2456  * \param theShapeType - type of sub-shape of theShape
2457  * \param theState - required state
2458  * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found sub-shapes
2459  */
2460 //=======================================================================
2461 Handle(TColStd_HSequenceOfInteger)
2462     GEOMImpl_IShapesOperations::GetShapesOnShapeIDs
2463                             (const Handle(GEOM_Object)& theCheckShape,
2464                              const Handle(GEOM_Object)& theShape,
2465                              const Standard_Integer theShapeType,
2466                              GEOMAlgo_State theState)
2467 {
2468   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
2469     getShapesOnShapeIDs (theCheckShape, theShape, theShapeType, theState);
2470
2471   if ( aSeqOfIDs.IsNull()  || aSeqOfIDs->Length() == 0 )
2472     return NULL;
2473
2474   // The GetShapesOnShape() doesn't change object so no new function is required.
2475   Handle(GEOM_Function) aFunction =
2476     GEOM::GetCreatedLast(theShape,theCheckShape)->GetLastFunction();
2477
2478   // Make a Python command
2479   GEOM::TPythonDump(aFunction)
2480     << "listShapesOnBoxIDs = geompy.GetShapesOnShapeIDs("
2481     << theCheckShape << ", "
2482     << theShape << ", "
2483     << TopAbs_ShapeEnum(theShapeType) << ", "
2484     << theState << ")";
2485
2486   SetErrorCode(OK);
2487   return aSeqOfIDs;
2488 }
2489
2490 //=======================================================================
2491 //function : GetShapesOnShape
2492 /*!
2493  * \brief Find sub-shapes complying with given status about surface
2494  * \param theCheckShape - the shape to check state of sub-shapes against
2495  * \param theShape - the shape to explore
2496  * \param theShapeType - type of sub-shape of theShape
2497  * \param theState - required state
2498  * \retval Handle(TColStd_HSequenceOfTransient) - found sub-shapes
2499  */
2500 //=======================================================================
2501 Handle(TColStd_HSequenceOfTransient)
2502   GEOMImpl_IShapesOperations::GetShapesOnShape
2503                              (const Handle(GEOM_Object)& theCheckShape,
2504                               const Handle(GEOM_Object)&  theShape,
2505                               const Standard_Integer theShapeType,
2506                               GEOMAlgo_State theState)
2507 {
2508   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
2509     getShapesOnShapeIDs (theCheckShape, theShape, theShapeType, theState);
2510   if ( aSeqOfIDs.IsNull()  || aSeqOfIDs->Length() == 0 )
2511     return NULL;
2512
2513   // Find objects by indices
2514   TCollection_AsciiString anAsciiList;
2515   Handle(TColStd_HSequenceOfTransient) aSeq;
2516   aSeq = getObjectsShapesOn( theShape, aSeqOfIDs, anAsciiList );
2517
2518   if ( aSeq.IsNull() || aSeq->IsEmpty() )
2519     return NULL;
2520
2521   // Make a Python command
2522
2523   Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast( aSeq->Value( 1 ));
2524   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
2525
2526   GEOM::TPythonDump(aFunction)
2527     << "[" << anAsciiList.ToCString() << "] = geompy.GetShapesOnShape("
2528     << theCheckShape << ", "
2529     << theShape << ", "
2530     << TopAbs_ShapeEnum(theShapeType) << ", "
2531     << theState << ")";
2532
2533   SetErrorCode(OK);
2534   return aSeq;
2535 }
2536
2537 //=======================================================================
2538 //function : GetShapesOnShapeAsCompound
2539 //=======================================================================
2540 Handle(GEOM_Object) GEOMImpl_IShapesOperations::GetShapesOnShapeAsCompound
2541                              (const Handle(GEOM_Object)& theCheckShape,
2542                               const Handle(GEOM_Object)&  theShape,
2543                               const Standard_Integer theShapeType,
2544                               GEOMAlgo_State theState)
2545 {
2546   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
2547     getShapesOnShapeIDs (theCheckShape, theShape, theShapeType, theState);
2548
2549   if ( aSeqOfIDs.IsNull()  || aSeqOfIDs->Length() == 0 )
2550     return NULL;
2551
2552   // Find objects by indices
2553   TCollection_AsciiString anAsciiList;
2554   Handle(TColStd_HSequenceOfTransient) aSeq;
2555   aSeq = getObjectsShapesOn( theShape, aSeqOfIDs, anAsciiList );
2556
2557   if ( aSeq.IsNull() || aSeq->IsEmpty() )
2558     return NULL;
2559
2560   TopoDS_Compound aCompound;
2561   BRep_Builder B;
2562   B.MakeCompound(aCompound);
2563   int i = 1;
2564   for(; i<=aSeq->Length(); i++) {
2565     Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast(aSeq->Value(i));
2566     TopoDS_Shape aShape_i = anObj->GetValue();
2567     B.Add(aCompound,aShape_i);
2568   }
2569
2570   //Add a new result object
2571   Handle(GEOM_Object) aRes = GetEngine()->AddObject(GetDocID(), GEOM_SHAPES_ON_SHAPE);
2572   Handle(GEOM_Function) aFunction =
2573     aRes->AddFunction(GEOMImpl_ShapeDriver::GetID(), SHAPES_ON_SHAPE);
2574   aFunction->SetValue(aCompound);
2575
2576   GEOM::TPythonDump(aFunction)
2577     << aRes << " = geompy.GetShapesOnShapeAsCompound("
2578     << theCheckShape << ", "
2579     << theShape << ", "
2580     << TopAbs_ShapeEnum(theShapeType) << ", "
2581     << theState << ")";
2582
2583   SetErrorCode(OK);
2584
2585   return aRes;
2586 }
2587
2588 //=======================================================================
2589 //function : getShapesOnSurfaceIDs
2590   /*!
2591    * \brief Find IDs of sub-shapes complying with given status about surface
2592     * \param theSurface - the surface to check state of sub-shapes against
2593     * \param theShape - the shape to explore
2594     * \param theShapeType - type of sub-shape of theShape
2595     * \param theState - required state
2596     * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found sub-shapes
2597    */
2598 //=======================================================================
2599 Handle(TColStd_HSequenceOfInteger)
2600   GEOMImpl_IShapesOperations::getShapesOnSurfaceIDs(const Handle(Geom_Surface)& theSurface,
2601                                                     const TopoDS_Shape&         theShape,
2602                                                     TopAbs_ShapeEnum            theShapeType,
2603                                                     GEOMAlgo_State              theState)
2604 {
2605   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs;
2606
2607   // Check presence of triangulation, build if need
2608   if (!CheckTriangulation(theShape)) {
2609     SetErrorCode("Cannot build triangulation on the shape");
2610     return aSeqOfIDs;
2611   }
2612
2613   // BEGIN: Mantis issue 0020961: Error on a pipe T-Shape
2614   // Compute tolerance
2615   Standard_Real T, VertMax = -RealLast();
2616   try {
2617 #if OCC_VERSION_LARGE > 0x06010000
2618     OCC_CATCH_SIGNALS;
2619 #endif
2620     for (TopExp_Explorer ExV (theShape, TopAbs_VERTEX); ExV.More(); ExV.Next()) {
2621       TopoDS_Vertex Vertex = TopoDS::Vertex(ExV.Current());
2622       T = BRep_Tool::Tolerance(Vertex);
2623       if (T > VertMax)
2624         VertMax = T;
2625     }
2626   }
2627   catch (Standard_Failure) {
2628     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2629     SetErrorCode(aFail->GetMessageString());
2630     return aSeqOfIDs;
2631   }
2632   // END: Mantis issue 0020961
2633
2634   // Call algo
2635   GEOMAlgo_FinderShapeOn1 aFinder;
2636   //Standard_Real aTol = 0.0001; // default value
2637   Standard_Real aTol = VertMax; // Mantis issue 0020961
2638
2639   aFinder.SetShape(theShape);
2640   aFinder.SetTolerance(aTol);
2641   aFinder.SetSurface(theSurface);
2642   aFinder.SetShapeType(theShapeType);
2643   aFinder.SetState(theState);
2644
2645   // Sets the minimal number of inner points for the faces that do not have own
2646   // inner points at all (for e.g. rectangular planar faces have just 2 triangles).
2647   // Default value=3
2648   aFinder.SetNbPntsMin(3);
2649   // Sets the maximal number of inner points for edges or faces.
2650   // It is usefull for the cases when this number is very big (e.g =2000) to improve
2651   // the performance. If this value =0, all inner points will be taken into account.
2652   // Default value=0
2653   aFinder.SetNbPntsMax(100);
2654
2655   aFinder.Perform();
2656
2657   // Interprete results
2658   Standard_Integer iErr = aFinder.ErrorStatus();
2659   // the detailed description of error codes is in GEOMAlgo_FinderShapeOn1.cxx
2660   if (iErr) {
2661     MESSAGE(" iErr : " << iErr);
2662     TCollection_AsciiString aMsg (" iErr : ");
2663     aMsg += TCollection_AsciiString(iErr);
2664     SetErrorCode(aMsg);
2665     return aSeqOfIDs;
2666   }
2667   Standard_Integer iWrn = aFinder.WarningStatus();
2668   // the detailed description of warning codes is in GEOMAlgo_FinderShapeOn1.cxx
2669   if (iWrn) {
2670     MESSAGE(" *** iWrn : " << iWrn);
2671   }
2672
2673   const TopTools_ListOfShape& listSS = aFinder.Shapes(); // the result
2674
2675   if (listSS.Extent() < 1) {
2676     //SetErrorCode("Not a single sub-shape of the requested type found on the given surface");
2677     SetErrorCode(NOT_FOUND_ANY); // NPAL18017
2678     return aSeqOfIDs;
2679   }
2680
2681   // Fill sequence of object IDs
2682   aSeqOfIDs = new TColStd_HSequenceOfInteger;
2683
2684   TopTools_IndexedMapOfShape anIndices;
2685   TopExp::MapShapes(theShape, anIndices);
2686
2687   TopTools_ListIteratorOfListOfShape itSub (listSS);
2688   for (int index = 1; itSub.More(); itSub.Next(), ++index) {
2689     int id = anIndices.FindIndex(itSub.Value());
2690     aSeqOfIDs->Append(id);
2691   }
2692
2693   return aSeqOfIDs;
2694 }
2695
2696 //=======================================================================
2697 //function : getObjectsShapesOn
2698 /*!
2699  * \brief Find shape objects and their entries by their ids
2700  * \param theShapeIDs - incoming shape ids
2701  * \param theShapeEntries - outgoing entries like "entry1, entry2, ..."
2702  * \retval Handle(TColStd_HSequenceOfTransient) - found shape objects
2703  */
2704 //=======================================================================
2705 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::
2706  getObjectsShapesOn(const Handle(GEOM_Object)&                theShape,
2707                     const Handle(TColStd_HSequenceOfInteger)& theShapeIDs,
2708                     TCollection_AsciiString &                 theShapeEntries)
2709 {
2710   Handle(TColStd_HSequenceOfTransient) aSeq;
2711
2712   if ( !theShapeIDs.IsNull() && theShapeIDs->Length() > 0 )
2713   {
2714     aSeq = new TColStd_HSequenceOfTransient;
2715     Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
2716     TCollection_AsciiString anEntry;
2717     for ( int i = 1; i <= theShapeIDs->Length(); ++i )
2718     {
2719       anArray->SetValue(1, theShapeIDs->Value( i ));
2720       Handle(GEOM_Object) anObj = GetEngine()->AddSubShape(theShape, anArray);
2721       aSeq->Append( anObj );
2722
2723       TDF_Tool::Entry(anObj->GetEntry(), anEntry);
2724       if ( i != 1 ) theShapeEntries += ",";
2725       theShapeEntries += anEntry;
2726     }
2727   }
2728   return aSeq;
2729 }
2730
2731 //=======================================================================
2732 //function : getShapesOnSurface
2733 /*!
2734    * \brief Find sub-shapes complying with given status about surface
2735     * \param theSurface - the surface to check state of sub-shapes against
2736     * \param theShape - the shape to explore
2737     * \param theShapeType - type of sub-shape of theShape
2738     * \param theState - required state
2739     * \param theShapeEntries - outgoing entries like "entry1, entry2, ..."
2740     * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found sub-shapes
2741  */
2742 //=======================================================================
2743 Handle(TColStd_HSequenceOfTransient)
2744     GEOMImpl_IShapesOperations::getShapesOnSurface(const Handle(Geom_Surface)& theSurface,
2745                                                    const Handle(GEOM_Object)&  theShape,
2746                                                    TopAbs_ShapeEnum            theShapeType,
2747                                                    GEOMAlgo_State              theState,
2748                                                    TCollection_AsciiString &   theShapeEntries)
2749 {
2750   // Find sub-shapes ids
2751   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
2752     getShapesOnSurfaceIDs (theSurface, theShape->GetValue(), theShapeType, theState);
2753   if ( aSeqOfIDs.IsNull()  || aSeqOfIDs->Length() == 0 )
2754     return NULL;
2755
2756   return getObjectsShapesOn( theShape, aSeqOfIDs, theShapeEntries );
2757 }
2758
2759 //=============================================================================
2760 /*!
2761  *  GetShapesOnPlane
2762  */
2763 //=============================================================================
2764 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetShapesOnPlane
2765                                         (const Handle(GEOM_Object)& theShape,
2766                                          const Standard_Integer     theShapeType,
2767                                          const Handle(GEOM_Object)& theAx1,
2768                                          const GEOMAlgo_State       theState)
2769 {
2770   SetErrorCode(KO);
2771
2772   if (theShape.IsNull() || theAx1.IsNull()) return NULL;
2773
2774   TopoDS_Shape aShape = theShape->GetValue();
2775   TopoDS_Shape anAx1  = theAx1->GetValue();
2776
2777   if (aShape.IsNull() || anAx1.IsNull()) return NULL;
2778
2779   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
2780   if ( !checkTypeShapesOn( theShapeType ))
2781     return NULL;
2782
2783   // Create plane
2784   Handle(Geom_Surface) aPlane = makePlane( anAx1 );
2785   if ( aPlane.IsNull() )
2786     return NULL;
2787
2788   // Find objects
2789   TCollection_AsciiString anAsciiList;
2790   Handle(TColStd_HSequenceOfTransient) aSeq;
2791   aSeq = getShapesOnSurface( aPlane, theShape, aShapeType, theState, anAsciiList );
2792   if ( aSeq.IsNull() || aSeq->Length() == 0 )
2793     return NULL;
2794
2795   // Make a Python command
2796
2797   Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast( aSeq->Value( 1 ));
2798   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
2799
2800   GEOM::TPythonDump(aFunction) << "[" << anAsciiList.ToCString()
2801     << "] = geompy.GetShapesOnPlane(" << theShape << ", "
2802       << aShapeType << ", " << theAx1 << ", " << theState << ")";
2803
2804   SetErrorCode(OK);
2805   return aSeq;
2806 }
2807
2808 //=============================================================================
2809 /*!
2810  *  GetShapesOnPlaneWithLocation
2811  */
2812 //=============================================================================
2813 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetShapesOnPlaneWithLocation
2814                                         (const Handle(GEOM_Object)& theShape,
2815                                          const Standard_Integer     theShapeType,
2816                                          const Handle(GEOM_Object)& theAx1,
2817                                          const Handle(GEOM_Object)& thePnt,
2818                                          const GEOMAlgo_State       theState)
2819 {
2820   SetErrorCode(KO);
2821
2822   if (theShape.IsNull() || theAx1.IsNull() || thePnt.IsNull()) return NULL;
2823
2824   TopoDS_Shape aShape = theShape->GetValue();
2825   TopoDS_Shape anAx1  = theAx1->GetValue();
2826   TopoDS_Shape anPnt = thePnt->GetValue();
2827
2828   if (aShape.IsNull() || anAx1.IsNull() || anPnt.IsNull()) return NULL;
2829
2830   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
2831   if ( !checkTypeShapesOn( theShapeType ))
2832     return NULL;
2833
2834   // Create plane
2835   if ( anAx1.ShapeType() != TopAbs_EDGE || anPnt.ShapeType() != TopAbs_VERTEX ) return NULL;
2836   TopoDS_Vertex V1, V2, V3;
2837   TopoDS_Edge anEdge = TopoDS::Edge(anAx1);
2838   TopExp::Vertices(anEdge, V1, V2, Standard_True);
2839
2840   if (V1.IsNull() || V2.IsNull()) {
2841     SetErrorCode("Bad edge given for the plane normal vector");
2842     return NULL;
2843   }
2844   V3 = TopoDS::Vertex(anPnt);
2845
2846   if(V3.IsNull()) {
2847     SetErrorCode("Bad vertex given for the plane location");
2848       return NULL;
2849   }
2850   gp_Pnt aLoc = BRep_Tool::Pnt(V3);
2851   gp_Vec aVec(BRep_Tool::Pnt(V1),BRep_Tool::Pnt(V2));
2852
2853   if (aVec.Magnitude() < Precision::Confusion()) {
2854     SetErrorCode("Vector with null magnitude given");
2855     return NULL;
2856   }
2857   Handle(Geom_Surface) aPlane = new Geom_Plane(aLoc, aVec);
2858
2859   if ( aPlane.IsNull() )
2860     return NULL;
2861
2862   // Find objects
2863   TCollection_AsciiString anAsciiList;
2864   Handle(TColStd_HSequenceOfTransient) aSeq;
2865   aSeq = getShapesOnSurface( aPlane, theShape, aShapeType, theState, anAsciiList );
2866   if ( aSeq.IsNull() || aSeq->Length() == 0 )
2867     return NULL;
2868
2869   // Make a Python command
2870
2871   Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast( aSeq->Value( 1 ));
2872   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
2873
2874   GEOM::TPythonDump(aFunction) << "[" << anAsciiList.ToCString()
2875     << "] = geompy.GetShapesOnPlaneWithLocation(" << theShape << ", "
2876     << aShapeType << ", " << theAx1 << ", "<< thePnt <<", " << theState << ")";
2877
2878   SetErrorCode(OK);
2879   return aSeq;
2880 }
2881
2882 //=============================================================================
2883 /*!
2884  *  GetShapesOnCylinder
2885  */
2886 //=============================================================================
2887 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetShapesOnCylinder
2888                                           (const Handle(GEOM_Object)& theShape,
2889                                            const Standard_Integer     theShapeType,
2890                                            const Handle(GEOM_Object)& theAxis,
2891                                            const Standard_Real        theRadius,
2892                                            const GEOMAlgo_State       theState)
2893 {
2894   SetErrorCode(KO);
2895
2896   if (theShape.IsNull() || theAxis.IsNull()) return NULL;
2897
2898   TopoDS_Shape aShape = theShape->GetValue();
2899   TopoDS_Shape anAxis = theAxis->GetValue();
2900
2901   if (aShape.IsNull() || anAxis.IsNull()) return NULL;
2902
2903   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
2904   if ( !checkTypeShapesOn( aShapeType ))
2905     return NULL;
2906
2907   // Create a cylinder surface
2908   Handle(Geom_Surface) aCylinder = makeCylinder( anAxis, theRadius );
2909   if ( aCylinder.IsNull() )
2910     return NULL;
2911
2912   // Find objects
2913   TCollection_AsciiString anAsciiList;
2914   Handle(TColStd_HSequenceOfTransient) aSeq;
2915   aSeq = getShapesOnSurface( aCylinder, theShape, aShapeType, theState, anAsciiList );
2916   if ( aSeq.IsNull() || aSeq->Length() == 0 )
2917     return NULL;
2918
2919   // Make a Python command
2920
2921   Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast( aSeq->Value( 1 ));
2922   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
2923
2924   GEOM::TPythonDump(aFunction) << "[" << anAsciiList.ToCString()
2925     << "] = geompy.GetShapesOnCylinder(" << theShape << ", " << aShapeType
2926       << ", " << theAxis << ", " << theRadius << ", " << theState << ")";
2927
2928   SetErrorCode(OK);
2929   return aSeq;
2930 }
2931
2932 //=============================================================================
2933 /*!
2934  *  GetShapesOnCylinderWithLocation
2935  */
2936 //=============================================================================
2937 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetShapesOnCylinderWithLocation
2938                                           (const Handle(GEOM_Object)& theShape,
2939                                            const Standard_Integer     theShapeType,
2940                                            const Handle(GEOM_Object)& theAxis,
2941                                            const Handle(GEOM_Object)& thePnt,
2942                                            const Standard_Real        theRadius,
2943                                            const GEOMAlgo_State       theState)
2944 {
2945   SetErrorCode(KO);
2946
2947   if (theShape.IsNull() || theAxis.IsNull() || thePnt.IsNull()) return NULL;
2948
2949   TopoDS_Shape aShape = theShape->GetValue();
2950   TopoDS_Shape anAxis = theAxis->GetValue();
2951   TopoDS_Shape aPnt   = thePnt->GetValue();
2952
2953   if (aShape.IsNull() || anAxis.IsNull() || aPnt.IsNull()) return NULL;
2954
2955   if (aPnt.ShapeType() != TopAbs_VERTEX )
2956   {
2957     SetErrorCode("Bottom location point must be vertex");
2958     return NULL;
2959   }
2960
2961   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
2962   if ( !checkTypeShapesOn( aShapeType ))
2963     return NULL;
2964
2965   // Create a cylinder surface
2966   Handle(Geom_Surface) aCylinder = makeCylinder( anAxis, theRadius );
2967   if ( aCylinder.IsNull() )
2968     return NULL;
2969
2970   // translate the surface
2971   Handle(Geom_CylindricalSurface) aCylSurface =
2972     Handle(Geom_CylindricalSurface)::DownCast( aCylinder );
2973   if ( aCylSurface.IsNull() )
2974   {
2975     SetErrorCode("Unexpected surface type instead of Geom_CylindricalSurface");
2976     return NULL;
2977   }
2978   gp_Pnt fromLoc = aCylSurface->Cylinder().Location();
2979   gp_Pnt toLoc   = BRep_Tool::Pnt( TopoDS::Vertex( aPnt ));
2980   aCylinder->Translate( fromLoc, toLoc );
2981
2982   // Find objects
2983   TCollection_AsciiString anAsciiList;
2984   Handle(TColStd_HSequenceOfTransient) aSeq;
2985   aSeq = getShapesOnSurface( aCylinder, theShape, aShapeType, theState, anAsciiList );
2986   if ( aSeq.IsNull() || aSeq->Length() == 0 )
2987     return NULL;
2988
2989   // Make a Python command
2990
2991   Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast( aSeq->Value( 1 ));
2992   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
2993
2994   GEOM::TPythonDump(aFunction)
2995     << "[" << anAsciiList.ToCString()
2996     << "] = geompy.GetShapesOnCylinderWithLocation(" << theShape << ", " << aShapeType << ", "
2997     << theAxis << ", " << thePnt << ", " << theRadius << ", " << theState << ")";
2998
2999   SetErrorCode(OK);
3000   return aSeq;
3001 }
3002
3003 //=============================================================================
3004 /*!
3005  *  GetShapesOnSphere
3006  */
3007 //=============================================================================
3008 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IShapesOperations::GetShapesOnSphere
3009                                           (const Handle(GEOM_Object)& theShape,
3010                                            const Standard_Integer     theShapeType,
3011                                            const Handle(GEOM_Object)& theCenter,
3012                                            const Standard_Real        theRadius,
3013                                            const GEOMAlgo_State       theState)
3014 {
3015   SetErrorCode(KO);
3016
3017   if (theShape.IsNull() || theCenter.IsNull()) return NULL;
3018
3019   TopoDS_Shape aShape  = theShape->GetValue();
3020   TopoDS_Shape aCenter = theCenter->GetValue();
3021
3022   if (aShape.IsNull() || aCenter.IsNull()) return NULL;
3023
3024   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
3025   if ( !checkTypeShapesOn( aShapeType ))
3026     return NULL;
3027
3028   // Center of the sphere
3029   if (aCenter.ShapeType() != TopAbs_VERTEX) return NULL;
3030   gp_Pnt aLoc = BRep_Tool::Pnt(TopoDS::Vertex(aCenter));
3031
3032   gp_Ax3 anAx3 (aLoc, gp::DZ());
3033   Handle(Geom_SphericalSurface) aSphere =
3034     new Geom_SphericalSurface(anAx3, theRadius);
3035
3036   // Find objects
3037   TCollection_AsciiString anAsciiList;
3038   Handle(TColStd_HSequenceOfTransient) aSeq;
3039   aSeq = getShapesOnSurface( aSphere, theShape, aShapeType, theState, anAsciiList );
3040   if ( aSeq.IsNull() || aSeq->Length() == 0 )
3041     return NULL;
3042
3043   // Make a Python command
3044
3045   Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast( aSeq->Value( 1 ));
3046   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
3047
3048   GEOM::TPythonDump(aFunction) << "[" << anAsciiList.ToCString()
3049     << "] = geompy.GetShapesOnSphere(" << theShape << ", " << aShapeType
3050       << ", " << theCenter << ", " << theRadius << ", " << theState << ")";
3051
3052   SetErrorCode(OK);
3053   return aSeq;
3054 }
3055
3056 //=============================================================================
3057 /*!
3058  *  GetShapesOnPlaneIDs
3059  */
3060 //=============================================================================
3061 Handle(TColStd_HSequenceOfInteger) GEOMImpl_IShapesOperations::GetShapesOnPlaneIDs
3062                                         (const Handle(GEOM_Object)& theShape,
3063                                          const Standard_Integer     theShapeType,
3064                                          const Handle(GEOM_Object)& theAx1,
3065                                          const GEOMAlgo_State       theState)
3066 {
3067   SetErrorCode(KO);
3068
3069   if (theShape.IsNull() || theAx1.IsNull()) return NULL;
3070
3071   TopoDS_Shape aShape = theShape->GetValue();
3072   TopoDS_Shape anAx1  = theAx1->GetValue();
3073
3074   if (aShape.IsNull() || anAx1.IsNull()) return NULL;
3075
3076   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
3077   if ( !checkTypeShapesOn( aShapeType ))
3078     return NULL;
3079
3080   // Create plane
3081   Handle(Geom_Surface) aPlane = makePlane( anAx1 );
3082   if ( aPlane.IsNull() )
3083     return NULL;
3084
3085   // Find object IDs
3086   Handle(TColStd_HSequenceOfInteger) aSeq;
3087   aSeq = getShapesOnSurfaceIDs( aPlane, aShape, aShapeType, theState );
3088
3089   // The GetShapesOnPlaneIDs() doesn't change object so no new function is required.
3090   Handle(GEOM_Function) aFunction = GEOM::GetCreatedLast(theShape,theAx1)->GetLastFunction();
3091
3092   // Make a Python command
3093   GEOM::TPythonDump(aFunction, /*append=*/true)
3094     << "listShapesOnPlane = geompy.GetShapesOnPlaneIDs"
3095     << "(" << theShape << "," << aShapeType << "," << theAx1 << "," << theState << ")";
3096
3097   SetErrorCode(OK);
3098   return aSeq;
3099 }
3100
3101 //=============================================================================
3102 /*!
3103  *  GetShapesOnPlaneWithLocationIDs
3104  */
3105 //=============================================================================
3106 Handle(TColStd_HSequenceOfInteger) GEOMImpl_IShapesOperations::GetShapesOnPlaneWithLocationIDs
3107                                         (const Handle(GEOM_Object)& theShape,
3108                                          const Standard_Integer     theShapeType,
3109                                          const Handle(GEOM_Object)& theAx1,
3110                                          const Handle(GEOM_Object)& thePnt,
3111                                          const GEOMAlgo_State       theState)
3112 {
3113   SetErrorCode(KO);
3114
3115   if (theShape.IsNull() || theAx1.IsNull() || thePnt.IsNull()) return NULL;
3116
3117   TopoDS_Shape aShape = theShape->GetValue();
3118   TopoDS_Shape anAx1  = theAx1->GetValue();
3119   TopoDS_Shape anPnt  = thePnt->GetValue();
3120
3121   if (aShape.IsNull() || anAx1.IsNull() || anPnt.IsNull()) return NULL;
3122
3123   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
3124   if ( !checkTypeShapesOn( aShapeType ))
3125     return NULL;
3126
3127   // Create plane
3128   if (anAx1.ShapeType() != TopAbs_EDGE || anPnt.ShapeType() != TopAbs_VERTEX) return NULL;
3129   TopoDS_Edge anEdge = TopoDS::Edge(anAx1);
3130   TopoDS_Vertex V1, V2, V3;
3131   TopExp::Vertices(anEdge, V1, V2, Standard_True);
3132   if (V1.IsNull() || V2.IsNull()) {
3133     SetErrorCode("Bad edge given for the plane normal vector");
3134     return NULL;
3135   }
3136   V3 = TopoDS::Vertex(anPnt);
3137   if(V3.IsNull()) {
3138     SetErrorCode("Bad vertex given for the plane location");
3139       return NULL;
3140   }
3141   gp_Pnt aLoc = BRep_Tool::Pnt(V3);
3142   gp_Vec aVec(BRep_Tool::Pnt(V1),BRep_Tool::Pnt(V2));
3143   if (aVec.Magnitude() < Precision::Confusion()) {
3144     SetErrorCode("Vector with null magnitude given");
3145     return NULL;
3146   }
3147
3148   Handle(Geom_Surface) aPlane = new Geom_Plane(aLoc, aVec);
3149   if ( aPlane.IsNull() )
3150     return NULL;
3151
3152   // Find object IDs
3153   Handle(TColStd_HSequenceOfInteger) aSeq;
3154   aSeq = getShapesOnSurfaceIDs( aPlane, aShape, aShapeType, theState );
3155
3156   // The GetShapesOnPlaneIDs() doesn't change object so no new function is required.
3157   Handle(GEOM_Function) aFunction = GEOM::GetCreatedLast(theShape,theAx1)->GetLastFunction();
3158
3159   // Make a Python command
3160   GEOM::TPythonDump(aFunction, /*append=*/true)
3161     << "listShapesOnPlane = geompy.GetShapesOnPlaneWithLocationIDs"
3162     << "(" << theShape << ", " << aShapeType << ", " << theAx1 << ", "<< thePnt << ", "  << theState << ")";
3163
3164   SetErrorCode(OK);
3165   return aSeq;
3166 }
3167
3168 //=============================================================================
3169 /*!
3170  *  GetShapesOnCylinderIDs
3171  */
3172 //=============================================================================
3173 Handle(TColStd_HSequenceOfInteger) GEOMImpl_IShapesOperations::GetShapesOnCylinderIDs
3174                                           (const Handle(GEOM_Object)& theShape,
3175                                            const Standard_Integer     theShapeType,
3176                                            const Handle(GEOM_Object)& theAxis,
3177                                            const Standard_Real        theRadius,
3178                                            const GEOMAlgo_State       theState)
3179 {
3180   SetErrorCode(KO);
3181
3182   if (theShape.IsNull() || theAxis.IsNull()) return NULL;
3183
3184   TopoDS_Shape aShape = theShape->GetValue();
3185   TopoDS_Shape anAxis = theAxis->GetValue();
3186
3187   if (aShape.IsNull() || anAxis.IsNull()) return NULL;
3188
3189   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
3190   if ( !checkTypeShapesOn( aShapeType ))
3191     return NULL;
3192
3193   // Create a cylinder surface
3194   Handle(Geom_Surface) aCylinder = makeCylinder( anAxis, theRadius );
3195   if ( aCylinder.IsNull() )
3196     return NULL;
3197
3198   // Find object IDs
3199   Handle(TColStd_HSequenceOfInteger) aSeq;
3200   aSeq = getShapesOnSurfaceIDs( aCylinder, aShape, aShapeType, theState );
3201
3202   // The GetShapesOnCylinder() doesn't change object so no new function is required.
3203   Handle(GEOM_Function) aFunction = GEOM::GetCreatedLast(theShape,theAxis)->GetLastFunction();
3204
3205   // Make a Python command
3206   GEOM::TPythonDump(aFunction, /*append=*/true)
3207     << "listShapesOnCylinder = geompy.GetShapesOnCylinderIDs"
3208     << "(" << theShape << ", " << aShapeType << ", " << theAxis << ", "
3209     << theRadius << ", " << theState << ")";
3210
3211   SetErrorCode(OK);
3212   return aSeq;
3213 }
3214
3215 //=============================================================================
3216 /*!
3217  *  GetShapesOnCylinderWithLocationIDs
3218  */
3219 //=============================================================================
3220 Handle(TColStd_HSequenceOfInteger) GEOMImpl_IShapesOperations::GetShapesOnCylinderWithLocationIDs
3221                                           (const Handle(GEOM_Object)& theShape,
3222                                            const Standard_Integer     theShapeType,
3223                                            const Handle(GEOM_Object)& theAxis,
3224                                            const Handle(GEOM_Object)& thePnt,
3225                                            const Standard_Real        theRadius,
3226                                            const GEOMAlgo_State       theState)
3227 {
3228   SetErrorCode(KO);
3229
3230   if (theShape.IsNull() || theAxis.IsNull() || thePnt.IsNull()) return NULL;
3231
3232   TopoDS_Shape aShape = theShape->GetValue();
3233   TopoDS_Shape anAxis = theAxis->GetValue();
3234   TopoDS_Shape aPnt   = thePnt->GetValue();
3235
3236   if (aShape.IsNull() || anAxis.IsNull() || aPnt.IsNull()) return NULL;
3237
3238   if (aPnt.ShapeType() != TopAbs_VERTEX )
3239   {
3240     SetErrorCode("Bottom location point must be vertex");
3241     return NULL;
3242   }
3243
3244   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
3245   if ( !checkTypeShapesOn( aShapeType ))
3246     return NULL;
3247
3248   // Create a cylinder surface
3249   Handle(Geom_Surface) aCylinder = makeCylinder( anAxis, theRadius );
3250   if ( aCylinder.IsNull() )
3251     return NULL;
3252
3253   // translate the surface
3254   Handle(Geom_CylindricalSurface) aCylSurface =
3255     Handle(Geom_CylindricalSurface)::DownCast( aCylinder );
3256   if ( aCylSurface.IsNull() )
3257   {
3258     SetErrorCode("Unexpected surface type instead of Geom_CylindricalSurface");
3259     return NULL;
3260   }
3261   gp_Pnt fromLoc = aCylSurface->Cylinder().Location();
3262   gp_Pnt toLoc   = BRep_Tool::Pnt( TopoDS::Vertex( aPnt ));
3263   aCylinder->Translate( fromLoc, toLoc );
3264
3265   // Find object IDs
3266   Handle(TColStd_HSequenceOfInteger) aSeq;
3267   aSeq = getShapesOnSurfaceIDs( aCylinder, aShape, aShapeType, theState );
3268
3269   // The GetShapesOnCylinder() doesn't change object so no new function is required.
3270   Handle(GEOM_Function) aFunction =
3271     GEOM::GetCreatedLast(theShape, GEOM::GetCreatedLast(thePnt,theAxis))->GetLastFunction();
3272
3273   // Make a Python command
3274   GEOM::TPythonDump(aFunction, /*append=*/true)
3275     << "listShapesOnCylinder = geompy.GetShapesOnCylinderWithLocationIDs"
3276     << "(" << theShape << ", " << aShapeType << ", " << theAxis << ", "
3277     << thePnt << ", " << theRadius << ", " << theState << ")";
3278
3279   SetErrorCode(OK);
3280   return aSeq;
3281 }
3282
3283 //=============================================================================
3284 /*!
3285  *  GetShapesOnSphereIDs
3286  */
3287 //=============================================================================
3288 Handle(TColStd_HSequenceOfInteger) GEOMImpl_IShapesOperations::GetShapesOnSphereIDs
3289                                           (const Handle(GEOM_Object)& theShape,
3290                                            const Standard_Integer     theShapeType,
3291                                            const Handle(GEOM_Object)& theCenter,
3292                                            const Standard_Real        theRadius,
3293                                            const GEOMAlgo_State       theState)
3294 {
3295   SetErrorCode(KO);
3296
3297   if (theShape.IsNull() || theCenter.IsNull()) return NULL;
3298
3299   TopoDS_Shape aShape  = theShape->GetValue();
3300   TopoDS_Shape aCenter = theCenter->GetValue();
3301
3302   if (aShape.IsNull() || aCenter.IsNull()) return NULL;
3303
3304   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
3305   if ( !checkTypeShapesOn( aShapeType ))
3306     return NULL;
3307
3308   // Center of the sphere
3309   if (aCenter.ShapeType() != TopAbs_VERTEX) return NULL;
3310   gp_Pnt aLoc = BRep_Tool::Pnt(TopoDS::Vertex(aCenter));
3311
3312   gp_Ax3 anAx3 (aLoc, gp::DZ());
3313   Handle(Geom_SphericalSurface) aSphere =
3314     new Geom_SphericalSurface(anAx3, theRadius);
3315
3316   // Find object IDs
3317   Handle(TColStd_HSequenceOfInteger) aSeq;
3318   aSeq = getShapesOnSurfaceIDs( aSphere, aShape, aShapeType, theState );
3319
3320   // The GetShapesOnSphere() doesn't change object so no new function is required.
3321   Handle(GEOM_Function) aFunction = GEOM::GetCreatedLast(theShape,theCenter)->GetLastFunction();
3322
3323   // Make a Python command
3324   GEOM::TPythonDump(aFunction, /*append=*/true)
3325     << "listShapesOnCylinder = geompy.GetShapesOnCylinderIDs"
3326     << "(" << theShape << ", " << aShapeType << ", " << theCenter << ", "
3327     << theRadius << ", " << theState << ")";
3328
3329   SetErrorCode(OK);
3330   return aSeq;
3331 }
3332
3333 //=======================================================================
3334 //function : getShapesOnQuadrangleIDs
3335   /*!
3336    * \brief Find IDs of sub-shapes complying with given status about quadrangle
3337     * \param theShape - the shape to explore
3338     * \param theShapeType - type of sub-shape of theShape
3339     * \param theTopLeftPoint - top left quadrangle corner
3340     * \param theTopRigthPoint - top right quadrangle corner
3341     * \param theBottomLeftPoint - bottom left quadrangle corner
3342     * \param theBottomRigthPoint - bottom right quadrangle corner
3343     * \param theState - required state
3344     * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found sub-shapes
3345    */
3346 //=======================================================================
3347 Handle(TColStd_HSequenceOfInteger)
3348   GEOMImpl_IShapesOperations::getShapesOnQuadrangleIDs (const Handle(GEOM_Object)& theShape,
3349                                                         const Standard_Integer     theShapeType,
3350                                                         const Handle(GEOM_Object)& theTopLeftPoint,
3351                                                         const Handle(GEOM_Object)& theTopRigthPoint,
3352                                                         const Handle(GEOM_Object)& theBottomLeftPoint,
3353                                                         const Handle(GEOM_Object)& theBottomRigthPoint,
3354                                                         const GEOMAlgo_State       theState)
3355 {
3356   SetErrorCode(KO);
3357
3358   if ( theShape.IsNull() ||
3359        theTopLeftPoint.IsNull() ||
3360        theTopRigthPoint.IsNull() ||
3361        theBottomLeftPoint.IsNull() ||
3362        theBottomRigthPoint.IsNull() )
3363     return NULL;
3364
3365   TopoDS_Shape aShape = theShape->GetValue();
3366   TopoDS_Shape aTL = theTopLeftPoint->GetValue();
3367   TopoDS_Shape aTR = theTopRigthPoint->GetValue();
3368   TopoDS_Shape aBL = theBottomLeftPoint->GetValue();
3369   TopoDS_Shape aBR = theBottomRigthPoint->GetValue();
3370
3371   if (aShape.IsNull() ||
3372       aTL.IsNull() ||
3373       aTR.IsNull() ||
3374       aBL.IsNull() ||
3375       aBR.IsNull() ||
3376       aTL.ShapeType() != TopAbs_VERTEX ||
3377       aTR.ShapeType() != TopAbs_VERTEX ||
3378       aBL.ShapeType() != TopAbs_VERTEX ||
3379       aBR.ShapeType() != TopAbs_VERTEX )
3380     return NULL;
3381
3382   TopAbs_ShapeEnum aShapeType = TopAbs_ShapeEnum(theShapeType);
3383   if ( !checkTypeShapesOn( aShapeType ))
3384     return NULL;
3385
3386   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs;
3387
3388   // Check presence of triangulation, build if need
3389   if (!CheckTriangulation(aShape)) {
3390     SetErrorCode("Cannot build triangulation on the shape");
3391     return aSeqOfIDs;
3392   }
3393
3394   // Call algo
3395   gp_Pnt aPntTL = BRep_Tool::Pnt(TopoDS::Vertex(aTL));
3396   gp_Pnt aPntTR = BRep_Tool::Pnt(TopoDS::Vertex(aTR));
3397   gp_Pnt aPntBL = BRep_Tool::Pnt(TopoDS::Vertex(aBL));
3398   gp_Pnt aPntBR = BRep_Tool::Pnt(TopoDS::Vertex(aBR));
3399
3400   GEOMAlgo_FinderShapeOnQuad aFinder( aPntTL, aPntTR, aPntBL, aPntBR );
3401   Standard_Real aTol = 0.0001; // default value
3402
3403   aFinder.SetShape(aShape);
3404   aFinder.SetTolerance(aTol);
3405   //aFinder.SetSurface(theSurface);
3406   aFinder.SetShapeType(aShapeType);
3407   aFinder.SetState(theState);
3408
3409   // Sets the minimal number of inner points for the faces that do not have own
3410   // inner points at all (for e.g. rectangular planar faces have just 2 triangles).
3411   // Default value=3
3412   aFinder.SetNbPntsMin(3);
3413   // Sets the maximal number of inner points for edges or faces.
3414   // It is usefull for the cases when this number is very big (e.g =2000) to improve
3415   // the performance. If this value =0, all inner points will be taken into account.
3416   // Default value=0
3417   aFinder.SetNbPntsMax(100);
3418
3419   aFinder.Perform();
3420
3421   // Interprete results
3422   Standard_Integer iErr = aFinder.ErrorStatus();
3423   // the detailed description of error codes is in GEOMAlgo_FinderShapeOn1.cxx
3424   if (iErr) {
3425     MESSAGE(" iErr : " << iErr);
3426     TCollection_AsciiString aMsg (" iErr : ");
3427     aMsg += TCollection_AsciiString(iErr);
3428     SetErrorCode(aMsg);
3429     return aSeqOfIDs;
3430   }
3431   Standard_Integer iWrn = aFinder.WarningStatus();
3432   // the detailed description of warning codes is in GEOMAlgo_FinderShapeOn1.cxx
3433   if (iWrn) {
3434     MESSAGE(" *** iWrn : " << iWrn);
3435   }
3436
3437   const TopTools_ListOfShape& listSS = aFinder.Shapes(); // the result
3438
3439   if (listSS.Extent() < 1) {
3440     //SetErrorCode("Not a single sub-shape of the requested type found on the given surface");
3441     SetErrorCode(NOT_FOUND_ANY); // NPAL18017
3442     return aSeqOfIDs;
3443   }
3444
3445   // Fill sequence of object IDs
3446   aSeqOfIDs = new TColStd_HSequenceOfInteger;
3447
3448   TopTools_IndexedMapOfShape anIndices;
3449   TopExp::MapShapes(aShape, anIndices);
3450
3451   TopTools_ListIteratorOfListOfShape itSub (listSS);
3452   for (int index = 1; itSub.More(); itSub.Next(), ++index) {
3453     int id = anIndices.FindIndex(itSub.Value());
3454     aSeqOfIDs->Append(id);
3455   }
3456   return aSeqOfIDs;
3457 }
3458
3459 //=======================================================================
3460 //function : GetShapesOnQuadrangle
3461   /*!
3462    * \brief Find sub-shapes complying with given status about quadrangle
3463     * \param theShape - the shape to explore
3464     * \param theShapeType - type of sub-shape of theShape
3465     * \param theTopLeftPoint - top left quadrangle corner
3466     * \param theTopRigthPoint - top right quadrangle corner
3467     * \param theBottomLeftPoint - bottom left quadrangle corner
3468     * \param theBottomRigthPoint - bottom right quadrangle corner
3469     * \param theState - required state
3470     * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found sub-shapes
3471    */
3472 //=======================================================================
3473 Handle(TColStd_HSequenceOfTransient)
3474     GEOMImpl_IShapesOperations::GetShapesOnQuadrangle (const Handle(GEOM_Object)& theShape,
3475                                                        const Standard_Integer     theShapeType,
3476                                                        const Handle(GEOM_Object)& theTopLeftPoint,
3477                                                        const Handle(GEOM_Object)& theTopRigthPoint,
3478                                                        const Handle(GEOM_Object)& theBottomLeftPoint,
3479                                                        const Handle(GEOM_Object)& theBottomRigthPoint,
3480                                                        const GEOMAlgo_State       theState)
3481 {
3482   // Find indices
3483   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
3484     getShapesOnQuadrangleIDs( theShape,
3485                               theShapeType,
3486                               theTopLeftPoint,
3487                               theTopRigthPoint,
3488                               theBottomLeftPoint,
3489                               theBottomRigthPoint,
3490                               theState);
3491   if ( aSeqOfIDs.IsNull() || aSeqOfIDs->IsEmpty() )
3492     return NULL;
3493
3494   // Find objects by indices
3495   TCollection_AsciiString anAsciiList;
3496   Handle(TColStd_HSequenceOfTransient) aSeq;
3497   aSeq = getObjectsShapesOn( theShape, aSeqOfIDs, anAsciiList );
3498   if ( aSeq.IsNull() || aSeq->IsEmpty() )
3499     return NULL;
3500
3501   // Make a Python command
3502
3503   Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast( aSeq->Value( 1 ));
3504   Handle(GEOM_Function) aFunction = anObj->GetLastFunction();
3505
3506   GEOM::TPythonDump(aFunction)
3507     << "[" << anAsciiList.ToCString() << "] = geompy.GetShapesOnQuadrangle("
3508     << theShape << ", "
3509     << TopAbs_ShapeEnum(theShapeType) << ", "
3510     << theTopLeftPoint << ", "
3511     << theTopRigthPoint << ", "
3512     << theBottomLeftPoint << ", "
3513     << theBottomRigthPoint << ", "
3514     << theState << ")";
3515
3516   SetErrorCode(OK);
3517   return aSeq;
3518 }
3519
3520 //=======================================================================
3521 //function : GetShapesOnQuadrangleIDs
3522   /*!
3523    * \brief Find IDs of sub-shapes complying with given status about quadrangle
3524     * \param theShape - the shape to explore
3525     * \param theShapeType - type of sub-shape of theShape
3526     * \param theTopLeftPoint - top left quadrangle corner
3527     * \param theTopRigthPoint - top right quadrangle corner
3528     * \param theBottomLeftPoint - bottom left quadrangle corner
3529     * \param theBottomRigthPoint - bottom right quadrangle corner
3530     * \param theState - required state
3531     * \retval Handle(TColStd_HSequenceOfInteger) - IDs of found sub-shapes
3532    */
3533 //=======================================================================
3534 Handle(TColStd_HSequenceOfInteger)
3535   GEOMImpl_IShapesOperations::GetShapesOnQuadrangleIDs (const Handle(GEOM_Object)& theShape,
3536                                                         const Standard_Integer     theShapeType,
3537                                                         const Handle(GEOM_Object)& theTopLeftPoint,
3538                                                         const Handle(GEOM_Object)& theTopRigthPoint,
3539                                                         const Handle(GEOM_Object)& theBottomLeftPoint,
3540                                                         const Handle(GEOM_Object)& theBottomRigthPoint,
3541                                                         const GEOMAlgo_State       theState)
3542 {
3543   // Find indices
3544   Handle(TColStd_HSequenceOfInteger) aSeqOfIDs =
3545     getShapesOnQuadrangleIDs( theShape,
3546                               theShapeType,
3547                               theTopLeftPoint,
3548                               theTopRigthPoint,
3549                               theBottomLeftPoint,
3550                               theBottomRigthPoint,
3551                               theState);
3552   if ( aSeqOfIDs.IsNull() || aSeqOfIDs->IsEmpty() )
3553     return NULL;
3554
3555   // Make a Python command
3556
3557   // The GetShapesOnCylinder() doesn't change object so no new function is required.
3558   Handle(GEOM_Object) lastObj = GEOM::GetCreatedLast(theShape,theTopLeftPoint);
3559   lastObj = GEOM::GetCreatedLast(lastObj,theTopRigthPoint);
3560   lastObj = GEOM::GetCreatedLast(lastObj,theBottomRigthPoint);
3561   lastObj = GEOM::GetCreatedLast(lastObj,theBottomLeftPoint);
3562   Handle(GEOM_Function) aFunction = lastObj->GetLastFunction();
3563
3564   GEOM::TPythonDump(aFunction, /*append=*/true)
3565     << "listShapesOnQuadrangle = geompy.GetShapesOnQuadrangleIDs("
3566     << theShape << ", "
3567     << TopAbs_ShapeEnum(theShapeType) << ", "
3568     << theTopLeftPoint << ", "
3569     << theTopRigthPoint << ", "
3570     << theBottomLeftPoint << ", "
3571     << theBottomRigthPoint << ", "
3572     << theState << ")";
3573
3574   SetErrorCode(OK);
3575   return aSeqOfIDs;
3576 }
3577
3578 //=============================================================================
3579 /*!
3580  *  GetInPlaceOfShape
3581  */
3582 //=============================================================================
3583 static bool GetInPlaceOfShape (const Handle(GEOM_Function)& theWhereFunction,
3584                                const TopTools_IndexedMapOfShape& theWhereIndices,
3585                                const TopoDS_Shape& theWhat,
3586                                TColStd_ListOfInteger& theModifiedList)
3587 {
3588   if (theWhereFunction.IsNull() || theWhat.IsNull()) return false;
3589
3590   if (theWhereIndices.Contains(theWhat)) {
3591     // entity was not changed by the operation
3592     Standard_Integer aWhatIndex = theWhereIndices.FindIndex(theWhat);
3593     theModifiedList.Append(aWhatIndex);
3594     return true;
3595   }
3596
3597   // try to find in history
3598   TDF_Label aHistoryLabel = theWhereFunction->GetHistoryEntry(Standard_False);
3599
3600   // search in history for all argument shapes
3601   Standard_Boolean isFound = Standard_False;
3602   Standard_Boolean isGood = Standard_False;
3603
3604   TDF_LabelSequence aLabelSeq;
3605   theWhereFunction->GetDependency(aLabelSeq);
3606   Standard_Integer nbArg = aLabelSeq.Length();
3607
3608   for (Standard_Integer iarg = 1; iarg <= nbArg && !isFound; iarg++) {
3609
3610     TDF_Label anArgumentRefLabel = aLabelSeq.Value(iarg);
3611
3612     Handle(GEOM_Object) anArgumentObject = GEOM_Object::GetReferencedObject(anArgumentRefLabel);
3613     TopoDS_Shape anArgumentShape = anArgumentObject->GetValue();
3614
3615     TopTools_IndexedMapOfShape anArgumentIndices;
3616     TopExp::MapShapes(anArgumentShape, anArgumentIndices);
3617
3618     if (anArgumentIndices.Contains(theWhat)) {
3619       isFound = Standard_True;
3620       Standard_Integer aWhatIndex = anArgumentIndices.FindIndex(theWhat);
3621
3622       // Find corresponding label in history
3623       TDF_Label anArgumentHistoryLabel =
3624         theWhereFunction->GetArgumentHistoryEntry(anArgumentRefLabel, Standard_False);
3625       if (anArgumentHistoryLabel.IsNull()) {
3626         // Lost History of operation argument. Possibly, all its entities was removed.
3627         isGood = Standard_True;
3628       }
3629       else {
3630         TDF_Label aWhatHistoryLabel = anArgumentHistoryLabel.FindChild(aWhatIndex, Standard_False);
3631
3632         if (aWhatHistoryLabel.IsNull()) {
3633           // Removed entity ? Compound ? Compsolid ? Shell ? Wire
3634           isGood = Standard_False;
3635         } else {
3636           Handle(TDataStd_IntegerArray) anIntegerArray;
3637           if (!aWhatHistoryLabel.FindAttribute(TDataStd_IntegerArray::GetID(), anIntegerArray)) {
3638             //Error: Empty modifications history for the sought shape.
3639             isGood = Standard_False;
3640           }
3641           else {
3642             isGood = Standard_True;
3643             Standard_Integer imod, aModifLen = anIntegerArray->Array()->Length();
3644             for (imod = 1; imod <= aModifLen; imod++) {
3645               theModifiedList.Append(anIntegerArray->Array()->Value(imod));
3646             }
3647           }
3648         }
3649       }
3650     }
3651   }
3652
3653   isFound = isGood;
3654
3655   if (!isFound) {
3656     // try compound/compsolid/shell/wire element by element
3657     bool isFoundAny = false;
3658     TopTools_MapOfShape mapShape;
3659
3660     if (theWhat.ShapeType() == TopAbs_COMPOUND ||
3661         theWhat.ShapeType() == TopAbs_COMPSOLID) {
3662       // recursive processing of compound/compsolid
3663       TopoDS_Iterator anIt (theWhat, Standard_True, Standard_True);
3664       for (; anIt.More(); anIt.Next()) {
3665         if (mapShape.Add(anIt.Value())) {
3666           TopoDS_Shape curWhat = anIt.Value();
3667           isFoundAny = GetInPlaceOfShape(theWhereFunction, theWhereIndices, curWhat, theModifiedList);
3668           if (isFoundAny) isFound = Standard_True;
3669         }
3670       }
3671     }
3672     else if (theWhat.ShapeType() == TopAbs_SHELL) {
3673       // try to replace a shell by its faces images
3674       TopExp_Explorer anExp (theWhat, TopAbs_FACE);
3675       for (; anExp.More(); anExp.Next()) {
3676         if (mapShape.Add(anExp.Current())) {
3677           TopoDS_Shape curWhat = anExp.Current();
3678           isFoundAny = GetInPlaceOfShape(theWhereFunction, theWhereIndices, curWhat, theModifiedList);
3679           if (isFoundAny) isFound = Standard_True;
3680         }
3681       }
3682     }
3683     else if (theWhat.ShapeType() == TopAbs_WIRE) {
3684       // try to replace a wire by its edges images
3685       TopExp_Explorer anExp (theWhat, TopAbs_EDGE);
3686       for (; anExp.More(); anExp.Next()) {
3687         if (mapShape.Add(anExp.Current())) {
3688           TopoDS_Shape curWhat = anExp.Current();
3689           isFoundAny = GetInPlaceOfShape(theWhereFunction, theWhereIndices, curWhat, theModifiedList);
3690           if (isFoundAny) isFound = Standard_True;
3691         }
3692       }
3693     }
3694     else {
3695       // Removed entity
3696     }
3697   }
3698
3699   return isFound;
3700 }
3701
3702 //=============================================================================
3703 /*!
3704  *  GetShapeProperties
3705  */
3706 //=============================================================================
3707 void GEOMImpl_IShapesOperations::GetShapeProperties( const TopoDS_Shape aShape, Standard_Real tab[],
3708                                                      gp_Pnt & aVertex )
3709 {
3710   GProp_GProps theProps;
3711   gp_Pnt aCenterMass;
3712   //TopoDS_Shape aPntShape;
3713   Standard_Real aShapeSize;
3714
3715   if    (aShape.ShapeType() == TopAbs_VERTEX) aCenterMass = BRep_Tool::Pnt( TopoDS::Vertex( aShape ) );
3716   else if (aShape.ShapeType() == TopAbs_EDGE) BRepGProp::LinearProperties(aShape,  theProps);
3717   else if (aShape.ShapeType() == TopAbs_FACE) BRepGProp::SurfaceProperties(aShape, theProps);
3718   else                                        BRepGProp::VolumeProperties(aShape,  theProps);
3719
3720   if (aShape.ShapeType() == TopAbs_VERTEX)
3721     aShapeSize = 1;
3722   else {
3723     aCenterMass = theProps.CentreOfMass();
3724     aShapeSize  = theProps.Mass();
3725   }
3726
3727 //   aPntShape = BRepBuilderAPI_MakeVertex(aCenterMass).Shape();
3728 //   aVertex   = BRep_Tool::Pnt( TopoDS::Vertex( aPntShape ) );
3729   aVertex = aCenterMass;
3730   tab[0] = aVertex.X();
3731   tab[1] = aVertex.Y();
3732   tab[2] = aVertex.Z();
3733   tab[3] = aShapeSize;
3734   return;
3735 }
3736
3737 namespace {
3738
3739   //================================================================================
3740   /*!
3741    * \brief Return normal to face at extrema point
3742    */
3743   //================================================================================
3744
3745   gp_Vec GetNormal (const TopoDS_Face& face, const BRepExtrema_DistShapeShape& extrema)
3746   {
3747     gp_Vec defaultNorm(1,0,0); // to have same normals on different faces
3748     try {
3749       // get UV at extrema point
3750       Standard_Real u,v, f,l;
3751       switch ( extrema.SupportTypeShape2(1) ) {
3752       case BRepExtrema_IsInFace: {
3753         extrema.ParOnFaceS2(1, u, v );
3754         break;
3755       }
3756       case BRepExtrema_IsOnEdge: {
3757         TopoDS_Edge edge = TopoDS::Edge( extrema.SupportOnShape2(1));
3758         Handle(Geom2d_Curve) pcurve = BRep_Tool::CurveOnSurface( edge, face, f,l );
3759         extrema.ParOnEdgeS2( 1, u );
3760         gp_Pnt2d uv = pcurve->Value( u );
3761         u = uv.Coord(1);
3762         v = uv.Coord(2);
3763         break;
3764       }
3765       case BRepExtrema_IsVertex: return defaultNorm;
3766       }
3767       // get derivatives
3768       BRepAdaptor_Surface surface( face, false );
3769       gp_Vec du, dv; gp_Pnt p;
3770       surface.D1( u, v, p, du, dv );
3771
3772       return du ^ dv;
3773
3774     } catch (Standard_Failure ) {
3775     }
3776     return defaultNorm;
3777   }
3778 }
3779
3780 //================================================================================
3781 /*!
3782  * \brief Return type of shape for explode. In case of compound it will be a type of sub-shape.
3783  */
3784 //================================================================================
3785 TopAbs_ShapeEnum GEOMImpl_IShapesOperations::GetTypeOfSimplePart (const TopoDS_Shape& theShape)
3786 {
3787   TopAbs_ShapeEnum aType = theShape.ShapeType();
3788   if      (aType == TopAbs_VERTEX)                             return TopAbs_VERTEX;
3789   else if (aType == TopAbs_EDGE  || aType == TopAbs_WIRE)      return TopAbs_EDGE;
3790   else if (aType == TopAbs_FACE  || aType == TopAbs_SHELL)     return TopAbs_FACE;
3791   else if (aType == TopAbs_SOLID || aType == TopAbs_COMPSOLID) return TopAbs_SOLID;
3792   else if (aType == TopAbs_COMPOUND) {
3793     // Only the iType of the first shape in the compound is taken into account
3794     TopoDS_Iterator It (theShape, Standard_False, Standard_False);
3795     if (It.More()) {
3796       return GetTypeOfSimplePart(It.Value());
3797     }
3798   }
3799   return TopAbs_SHAPE;
3800 }
3801
3802 //=============================================================================
3803 /*!
3804  *  case GetInPlace:
3805  *  default:
3806  */
3807 //=============================================================================
3808 Handle(GEOM_Object) GEOMImpl_IShapesOperations::GetInPlace (Handle(GEOM_Object) theShapeWhere,
3809                                                             Handle(GEOM_Object) theShapeWhat)
3810 {
3811   SetErrorCode(KO);
3812
3813   if (theShapeWhere.IsNull() || theShapeWhat.IsNull()) return NULL;
3814
3815   TopoDS_Shape aWhere = theShapeWhere->GetValue();
3816   TopoDS_Shape aWhat  = theShapeWhat->GetValue();
3817   TopoDS_Shape aPntShape;
3818   TopoDS_Vertex aVertex;
3819
3820   if (aWhere.IsNull() || aWhat.IsNull()) {
3821     SetErrorCode("Error: aWhere and aWhat TopoDS_Shape are Null.");
3822     return NULL;
3823   }
3824
3825   Handle(GEOM_Function) aWhereFunction = theShapeWhere->GetLastFunction();
3826   if (aWhereFunction.IsNull()) {
3827     SetErrorCode("Error: aWhereFunction is Null.");
3828     return NULL;
3829   }
3830
3831   TopTools_IndexedMapOfShape aWhereIndices;
3832   TopExp::MapShapes(aWhere, aWhereIndices);
3833
3834   TopAbs_ShapeEnum iType = TopAbs_SOLID;
3835   Standard_Real    dl_l = 1e-3;
3836   Standard_Real    min_l, Tol_0D, Tol_1D, Tol_2D, Tol_3D, Tol_Mass;
3837   Standard_Real    aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
3838   Bnd_Box          BoundingBox;
3839   gp_Pnt           aPnt, aPnt_aWhat, tab_Pnt[2];
3840   GProp_GProps     aProps;
3841
3842   // Find the iType of the aWhat shape
3843   iType = GetTypeOfSimplePart(aWhat);
3844   if (iType == TopAbs_SHAPE) {
3845     SetErrorCode("Error: An attempt to extract a shape of not supported type.");
3846     return NULL;
3847   }
3848
3849   TopExp_Explorer Exp_aWhat  ( aWhat,  iType );
3850   TopExp_Explorer Exp_aWhere ( aWhere, iType );
3851   TopExp_Explorer Exp_Edge   ( aWhere, TopAbs_EDGE );
3852
3853   // Find the shortest edge in theShapeWhere shape
3854   BRepBndLib::Add(aWhere, BoundingBox);
3855   BoundingBox.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
3856   min_l = fabs(aXmax - aXmin);
3857   if( min_l < fabs(aYmax - aYmin) ) min_l = fabs(aYmax - aYmin);
3858   if( min_l < fabs(aZmax - aZmin) ) min_l = fabs(aZmax - aZmin);
3859   min_l /= dl_l;
3860   // Mantis issue 0020908 BEGIN
3861   if (!Exp_Edge.More()) {
3862     min_l = Precision::Confusion();
3863   }
3864   // Mantis issue 0020908 END
3865   for ( Standard_Integer nbEdge = 0; Exp_Edge.More(); Exp_Edge.Next(), nbEdge++ ) {
3866     TopExp_Explorer Exp_Vertex( Exp_Edge.Current(), TopAbs_VERTEX);
3867     for ( Standard_Integer nbVertex = 0; Exp_Vertex.More(); Exp_Vertex.Next(), nbVertex++ ) {
3868       aPnt = BRep_Tool::Pnt( TopoDS::Vertex( Exp_Vertex.Current() ) );
3869       tab_Pnt[nbVertex] = aPnt;
3870     }
3871     if ( ! tab_Pnt[0].IsEqual(tab_Pnt[1], dl_l) ) {
3872       BRepGProp::LinearProperties(Exp_Edge.Current(), aProps);
3873       if ( aProps.Mass() < min_l ) min_l = aProps.Mass();
3874     }
3875   }
3876
3877   // Compute tolerances
3878   Tol_0D = dl_l;
3879   Tol_1D = dl_l * min_l;
3880   Tol_2D = dl_l * ( min_l * min_l) * ( 2. + dl_l);
3881   Tol_3D = dl_l * ( min_l * min_l * min_l ) * ( 3. + (3 * dl_l) + (dl_l * dl_l) );
3882
3883   if (Tol_0D < Precision::Confusion()) Tol_0D = Precision::Confusion();
3884   if (Tol_1D < Precision::Confusion()) Tol_1D = Precision::Confusion();
3885   if (Tol_2D < Precision::Confusion()) Tol_2D = Precision::Confusion();
3886   if (Tol_3D < Precision::Confusion()) Tol_3D = Precision::Confusion();
3887
3888   Tol_Mass = Tol_3D;
3889   if ( iType == TopAbs_VERTEX )    Tol_Mass = Tol_0D;
3890   else if ( iType == TopAbs_EDGE ) Tol_Mass = Tol_1D;
3891   else if ( iType == TopAbs_FACE ) Tol_Mass = Tol_2D;
3892
3893   // Searching for the sub-shapes inside the ShapeWhere shape
3894   GEOMAlgo_GetInPlace aGIP;
3895   aGIP.SetTolerance(Tol_1D);
3896   aGIP.SetTolMass(Tol_Mass);
3897   aGIP.SetTolCG(Tol_1D);
3898
3899   aGIP.SetArgument(aWhat);
3900   aGIP.SetShapeWhere(aWhere);
3901
3902   aGIP.Perform();
3903   int iErr = aGIP.ErrorStatus();
3904   if (iErr) {
3905     SetErrorCode("Error in GEOMAlgo_GetInPlace");
3906     return NULL;
3907   }
3908
3909   // aGIP.IsFound() returns true only when the whole theShapeWhat
3910   // is found (as one shape or several parts). But we are also interested
3911   // in the partial result, that is why this check is commented.
3912   //if (!aGIP.IsFound()) {
3913   //  SetErrorCode(NOT_FOUND_ANY);
3914   //  return NULL;
3915   //}
3916
3917   const TopTools_DataMapOfShapeListOfShape& aDMSLS = aGIP.Images();
3918   if (!aDMSLS.IsBound(aWhat)) {
3919     SetErrorCode(NOT_FOUND_ANY);
3920     return NULL;
3921   }
3922
3923   // the list of shapes aLSA contains the shapes 
3924   // of the Shape For Search that corresponds 
3925   // to the Argument aWhat
3926   const TopTools_ListOfShape& aLSA = aDMSLS.Find(aWhat);
3927   if (aLSA.Extent() == 0) {
3928     SetErrorCode(NOT_FOUND_ANY); // Not found any Results
3929     return NULL;
3930   }
3931
3932   Handle(TColStd_HArray1OfInteger) aModifiedArray = new TColStd_HArray1OfInteger (1, aLSA.Extent());
3933   TopTools_ListIteratorOfListOfShape anIterModif (aLSA);
3934   for (Standard_Integer imod = 1; anIterModif.More(); anIterModif.Next(), imod++) {
3935     if (aWhereIndices.Contains(anIterModif.Value())) {
3936       aModifiedArray->SetValue(imod, aWhereIndices.FindIndex(anIterModif.Value()));
3937     }
3938     else {
3939       SetErrorCode("Error: wrong sub-shape returned");
3940       return NULL;
3941     }
3942   }
3943
3944   //Add a new object
3945   Handle(GEOM_Object) aResult = GetEngine()->AddSubShape(theShapeWhere, aModifiedArray);
3946   if (aResult.IsNull()) {
3947     SetErrorCode("Error in algorithm: result found, but cannot be returned.");
3948     return NULL;
3949   }
3950
3951   if (aModifiedArray->Length() > 1 || theShapeWhat->GetType() == GEOM_GROUP) {
3952     //Set a GROUP type
3953     aResult->SetType(GEOM_GROUP);
3954
3955     //Set a sub-shape type
3956     TopoDS_Shape aFirstFound = aLSA.First();
3957     TopAbs_ShapeEnum aShapeType = aFirstFound.ShapeType();
3958
3959     TDF_Label aFreeLabel = aResult->GetFreeLabel();
3960     TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)aShapeType);
3961   }
3962
3963   //Make a Python command
3964   Handle(GEOM_Function) aFunction = aResult->GetFunction(1);
3965
3966   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetInPlace("
3967     << theShapeWhere << ", " << theShapeWhat << ", True)";
3968
3969   SetErrorCode(OK);
3970   return aResult;
3971 }
3972
3973 //=============================================================================
3974 /*!
3975  *  case GetInPlaceOld:
3976  *  default:
3977  */
3978 //=============================================================================
3979 Handle(GEOM_Object) GEOMImpl_IShapesOperations::GetInPlaceOld (Handle(GEOM_Object) theShapeWhere,
3980                                                                Handle(GEOM_Object) theShapeWhat)
3981 {
3982   SetErrorCode(KO);
3983
3984   if (theShapeWhere.IsNull() || theShapeWhat.IsNull()) return NULL;
3985
3986   TopoDS_Shape aWhere = theShapeWhere->GetValue();
3987   TopoDS_Shape aWhat  = theShapeWhat->GetValue();
3988   TopoDS_Shape aPntShape;
3989   TopoDS_Vertex aVertex;
3990
3991   if (aWhere.IsNull() || aWhat.IsNull()) {
3992     SetErrorCode("Error: aWhere and aWhat TopoDS_Shape are Null.");
3993     return NULL;
3994   }
3995
3996   Handle(GEOM_Function) aWhereFunction = theShapeWhere->GetLastFunction();
3997   if (aWhereFunction.IsNull()) {
3998     SetErrorCode("Error: aWhereFunction is Null.");
3999     return NULL;
4000   }
4001
4002   TopTools_IndexedMapOfShape aWhereIndices;
4003   TopExp::MapShapes(aWhere, aWhereIndices);
4004
4005   TColStd_ListOfInteger aModifiedList;
4006   Standard_Integer aWhereIndex;
4007   Handle(TColStd_HArray1OfInteger) aModifiedArray;
4008   Handle(GEOM_Object) aResult;
4009
4010   bool isFound = false;
4011   TopAbs_ShapeEnum iType = TopAbs_SOLID;
4012   //Standard_Real    aWhat_Mass = 0., aWhere_Mass = 0.;
4013   Standard_Real    tab_aWhat[4],    tab_aWhere[4];
4014   Standard_Real    dl_l = 1e-3;
4015   Standard_Real    min_l, Tol_0D, Tol_1D, Tol_2D, Tol_3D, Tol_Mass;
4016   Standard_Real    aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
4017   Bnd_Box          BoundingBox;
4018   gp_Pnt           aPnt, aPnt_aWhat, tab_Pnt[2];
4019   GProp_GProps     aProps;
4020
4021   // Find the iType of the aWhat shape
4022   /*
4023   if      ( aWhat.ShapeType() == TopAbs_VERTEX )                                         iType = TopAbs_VERTEX;
4024   else if ( aWhat.ShapeType() == TopAbs_EDGE  || aWhat.ShapeType() == TopAbs_WIRE )      iType = TopAbs_EDGE;
4025   else if ( aWhat.ShapeType() == TopAbs_FACE  || aWhat.ShapeType() == TopAbs_SHELL )     iType = TopAbs_FACE;
4026   else if ( aWhat.ShapeType() == TopAbs_SOLID || aWhat.ShapeType() == TopAbs_COMPSOLID ) iType = TopAbs_SOLID;
4027   else if ( aWhat.ShapeType() == TopAbs_COMPOUND ) {
4028     // Only the iType of the first shape in the compound is taken into account
4029     TopoDS_Iterator It (aWhat, Standard_False, Standard_False);
4030     if ( !It.More() ) {
4031       SetErrorCode("Error: theShapeWhat is an empty COMPOUND.");
4032       return NULL;
4033     }
4034     TopAbs_ShapeEnum compType = It.Value().ShapeType();
4035     if      ( compType == TopAbs_VERTEX )                               iType = TopAbs_VERTEX;
4036     else if ( compType == TopAbs_EDGE  || compType == TopAbs_WIRE )     iType = TopAbs_EDGE;
4037     else if ( compType == TopAbs_FACE  || compType == TopAbs_SHELL)     iType = TopAbs_FACE;
4038     else if ( compType == TopAbs_SOLID || compType == TopAbs_COMPSOLID) iType = TopAbs_SOLID;
4039   }
4040   else {
4041     SetErrorCode("Error: An attempt to extract a shape of not supported type.");
4042     return NULL;
4043   }
4044   */
4045   iType = GetTypeOfSimplePart(aWhat);
4046   if (iType == TopAbs_SHAPE) {
4047     SetErrorCode("Error: An attempt to extract a shape of not supported type.");
4048     return NULL;
4049   }
4050
4051   TopExp_Explorer Exp_aWhat  ( aWhat,  iType );
4052   TopExp_Explorer Exp_aWhere ( aWhere, iType );
4053   TopExp_Explorer Exp_Edge   ( aWhere, TopAbs_EDGE );
4054
4055   // Find the shortest edge in theShapeWhere shape
4056   BRepBndLib::Add(aWhere, BoundingBox);
4057   BoundingBox.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
4058   min_l = fabs(aXmax - aXmin);
4059   if( min_l < fabs(aYmax - aYmin) ) min_l = fabs(aYmax - aYmin);
4060   if( min_l < fabs(aZmax - aZmin) ) min_l = fabs(aZmax - aZmin);
4061   min_l /= dl_l;
4062   // Mantis issue 0020908 BEGIN
4063   if (!Exp_Edge.More()) {
4064     min_l = Precision::Confusion();
4065   }
4066   // Mantis issue 0020908 END
4067   for ( Standard_Integer nbEdge = 0; Exp_Edge.More(); Exp_Edge.Next(), nbEdge++ ) {
4068     TopExp_Explorer Exp_Vertex( Exp_Edge.Current(), TopAbs_VERTEX);
4069     for ( Standard_Integer nbVertex = 0; Exp_Vertex.More(); Exp_Vertex.Next(), nbVertex++ ) {
4070       aPnt = BRep_Tool::Pnt( TopoDS::Vertex( Exp_Vertex.Current() ) );
4071       tab_Pnt[nbVertex] = aPnt;
4072     }
4073     if ( ! tab_Pnt[0].IsEqual(tab_Pnt[1], dl_l) ) {
4074       BRepGProp::LinearProperties(Exp_Edge.Current(), aProps);
4075       if ( aProps.Mass() < min_l ) min_l = aProps.Mass();
4076     }
4077   }
4078
4079   // Compute tolerances
4080   Tol_0D = dl_l;
4081   Tol_1D = dl_l * min_l;
4082   Tol_2D = dl_l * ( min_l * min_l) * ( 2. + dl_l);
4083   Tol_3D = dl_l * ( min_l * min_l * min_l ) * ( 3. + (3 * dl_l) + (dl_l * dl_l) );
4084
4085   if (Tol_0D < Precision::Confusion()) Tol_0D = Precision::Confusion();
4086   if (Tol_1D < Precision::Confusion()) Tol_1D = Precision::Confusion();
4087   if (Tol_2D < Precision::Confusion()) Tol_2D = Precision::Confusion();
4088   if (Tol_3D < Precision::Confusion()) Tol_3D = Precision::Confusion();
4089
4090   //if (Tol_1D > 1.0) Tol_1D = 1.0;
4091   //if (Tol_2D > 1.0) Tol_2D = 1.0;
4092   //if (Tol_3D > 1.0) Tol_3D = 1.0;
4093
4094   Tol_Mass = Tol_3D;
4095   if ( iType == TopAbs_VERTEX )    Tol_Mass = Tol_0D;
4096   else if ( iType == TopAbs_EDGE ) Tol_Mass = Tol_1D;
4097   else if ( iType == TopAbs_FACE ) Tol_Mass = Tol_2D;
4098
4099   // Compute the ShapeWhat Mass
4100   /*
4101   for ( ; Exp_aWhat.More(); Exp_aWhat.Next() ) {
4102     if ( iType == TopAbs_VERTEX ) {
4103       aWhat_Mass += 1;
4104       continue;
4105     }
4106     else if ( iType == TopAbs_EDGE ) BRepGProp::LinearProperties(Exp_aWhat.Current(),  aProps);
4107     else if ( iType == TopAbs_FACE ) BRepGProp::SurfaceProperties(Exp_aWhat.Current(), aProps);
4108     else                             BRepGProp::VolumeProperties(Exp_aWhat.Current(),  aProps);
4109     aWhat_Mass += aProps.Mass();
4110   }
4111   */
4112
4113   // Searching for the sub-shapes inside the ShapeWhere shape
4114   TopTools_MapOfShape map_aWhere;
4115   for ( Exp_aWhere.ReInit(); Exp_aWhere.More(); Exp_aWhere.Next() ) {
4116     if (!map_aWhere.Add(Exp_aWhere.Current()))
4117       continue; // skip repeated shape to avoid mass addition
4118     GetShapeProperties( Exp_aWhere.Current(), tab_aWhere, aPnt );
4119     for ( Exp_aWhat.ReInit(); Exp_aWhat.More(); Exp_aWhat.Next() ) {
4120       GetShapeProperties( Exp_aWhat.Current(), tab_aWhat, aPnt_aWhat );
4121       if ( fabs(tab_aWhat[3] - tab_aWhere[3]) <= Tol_Mass && aPnt_aWhat.Distance(aPnt) <= Tol_1D )
4122         isFound = true;
4123       else {
4124         if ( (tab_aWhat[3] - tab_aWhere[3]) > Tol_Mass ) {
4125           aPntShape = BRepBuilderAPI_MakeVertex( aPnt ).Shape();
4126           aVertex   = TopoDS::Vertex( aPntShape );
4127           BRepExtrema_DistShapeShape aWhereDistance ( aVertex, Exp_aWhere.Current() );
4128           BRepExtrema_DistShapeShape aWhatDistance  ( aVertex, Exp_aWhat.Current() );
4129           if ( aWhereDistance.IsDone() && aWhatDistance.IsDone() &&
4130                fabs(aWhereDistance.Value() - aWhatDistance.Value()) <= Tol_1D )
4131           {
4132             // 0020162: "EDF 961 GEOM : Getinplace is getting additionnal orthogonal faces"
4133             // aVertex must be projected to the same point on Where and on What
4134             gp_Pnt pOnWhat  = aWhatDistance.PointOnShape2(1);
4135             gp_Pnt pOnWhere = aWhereDistance.PointOnShape2(1);
4136             isFound = ( pOnWhat.Distance(pOnWhere) <= Tol_1D );
4137             if ( isFound && iType == TopAbs_FACE )
4138             {
4139               // check normals at pOnWhat and pOnWhere
4140               const double angleTol = M_PI/180.;
4141               gp_Vec normToWhat  = GetNormal( TopoDS::Face(Exp_aWhat.Current()), aWhatDistance);
4142               gp_Vec normToWhere = GetNormal( TopoDS::Face(Exp_aWhere.Current()), aWhereDistance);
4143               if ( normToWhat * normToWhere < 0 )
4144                 normToWhat.Reverse();
4145               isFound = ( normToWhat.Angle( normToWhere ) < angleTol );
4146             }
4147           }
4148         }
4149       }
4150       if ( isFound ) {
4151         aWhereIndex = aWhereIndices.FindIndex(Exp_aWhere.Current());
4152         aModifiedList.Append(aWhereIndex);
4153         //aWhere_Mass += tab_aWhere[3];
4154         isFound = false;
4155         break;
4156       }
4157     }
4158     //if ( fabs( aWhat_Mass - aWhere_Mass ) <= Tol_Mass )
4159       //break;
4160   }
4161
4162   if (aModifiedList.Extent() == 0) { // Not found any Results
4163     SetErrorCode(NOT_FOUND_ANY);
4164     return NULL;
4165   }
4166
4167   aModifiedArray = new TColStd_HArray1OfInteger (1, aModifiedList.Extent());
4168   TColStd_ListIteratorOfListOfInteger anIterModif (aModifiedList);
4169   for (Standard_Integer imod = 1; anIterModif.More(); anIterModif.Next(), imod++)
4170     aModifiedArray->SetValue(imod, anIterModif.Value());
4171
4172   //Add a new object
4173   aResult = GetEngine()->AddSubShape(theShapeWhere, aModifiedArray);
4174   if (aResult.IsNull()) {
4175     SetErrorCode("Error in algorithm: result found, but cannot be returned.");
4176     return NULL;
4177   }
4178
4179   if (aModifiedArray->Length() > 1 || theShapeWhat->GetType() == GEOM_GROUP) {
4180     //Set a GROUP type
4181     aResult->SetType(GEOM_GROUP);
4182
4183     //Set a sub-shape type
4184     TopoDS_Shape aFirstFound = aWhereIndices.FindKey(aModifiedArray->Value(1));
4185     TopAbs_ShapeEnum aShapeType = aFirstFound.ShapeType();
4186
4187     TDF_Label aFreeLabel = aResult->GetFreeLabel();
4188     TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)aShapeType);
4189   }
4190
4191   //Make a Python command
4192   Handle(GEOM_Function) aFunction = aResult->GetFunction(1);
4193
4194   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetInPlace("
4195     << theShapeWhere << ", " << theShapeWhat << ", False)";
4196
4197   SetErrorCode(OK);
4198   return aResult;
4199 }
4200
4201 //=======================================================================
4202 //function : GetInPlaceByHistory
4203 //purpose  :
4204 //=======================================================================
4205 Handle(GEOM_Object) GEOMImpl_IShapesOperations::GetInPlaceByHistory
4206                                           (Handle(GEOM_Object) theShapeWhere,
4207                                            Handle(GEOM_Object) theShapeWhat)
4208 {
4209   SetErrorCode(KO);
4210
4211   if (theShapeWhere.IsNull() || theShapeWhat.IsNull()) return NULL;
4212
4213   TopoDS_Shape aWhere = theShapeWhere->GetValue();
4214   TopoDS_Shape aWhat  = theShapeWhat->GetValue();
4215
4216   if (aWhere.IsNull() || aWhat.IsNull()) return NULL;
4217
4218   Handle(GEOM_Function) aWhereFunction = theShapeWhere->GetLastFunction();
4219   if (aWhereFunction.IsNull()) return NULL;
4220
4221   //Fill array of indices
4222   TopTools_IndexedMapOfShape aWhereIndices;
4223   TopExp::MapShapes(aWhere, aWhereIndices);
4224
4225   // process shape
4226   TColStd_ListOfInteger aModifiedList;
4227   bool isFound = GetInPlaceOfShape(aWhereFunction, aWhereIndices, aWhat, aModifiedList);
4228
4229   if (!isFound || aModifiedList.Extent() < 1) {
4230     SetErrorCode("Error: No history found for the sought shape or its sub-shapes.");
4231     return NULL;
4232   }
4233
4234   Handle(TColStd_HArray1OfInteger) aModifiedArray =
4235     new TColStd_HArray1OfInteger (1, aModifiedList.Extent());
4236   TColStd_ListIteratorOfListOfInteger anIterModif (aModifiedList);
4237   for (Standard_Integer imod = 1; anIterModif.More(); anIterModif.Next(), imod++) {
4238     aModifiedArray->SetValue(imod, anIterModif.Value());
4239   }
4240
4241   //Add a new object
4242   Handle(GEOM_Object) aResult = GetEngine()->AddSubShape(theShapeWhere, aModifiedArray);
4243   if (aResult.IsNull()) {
4244     SetErrorCode("Error in algorithm: result found, but cannot be returned.");
4245     return NULL;
4246   }
4247
4248   if (aModifiedArray->Length() > 1) {
4249     //Set a GROUP type
4250     aResult->SetType(GEOM_GROUP);
4251
4252     //Set a sub-shape type
4253     TopoDS_Shape aFirstFound = aWhereIndices.FindKey(aModifiedArray->Value(1));
4254     TopAbs_ShapeEnum aShapeType = aFirstFound.ShapeType();
4255
4256     TDF_Label aFreeLabel = aResult->GetFreeLabel();
4257     TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)aShapeType);
4258   }
4259
4260   //Make a Python command
4261   Handle(GEOM_Function) aFunction = aResult->GetFunction(1);
4262
4263   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetInPlaceByHistory("
4264     << theShapeWhere << ", " << theShapeWhat << ")";
4265
4266   SetErrorCode(OK);
4267   return aResult;
4268 }
4269
4270 //=======================================================================
4271 //function : ShapeToDouble
4272 //purpose  : used by CompareShapes::operator()
4273 //=======================================================================
4274 std::pair<double, double> ShapeToDouble (const TopoDS_Shape& S, bool isOldSorting)
4275 {
4276   // Computing of CentreOfMass
4277   gp_Pnt GPoint;
4278   double Len;
4279
4280   if (S.ShapeType() == TopAbs_VERTEX) {
4281     GPoint = BRep_Tool::Pnt(TopoDS::Vertex(S));
4282     Len = (double)S.Orientation();
4283   }
4284   else {
4285     GProp_GProps GPr;
4286     // BEGIN: fix for Mantis issue 0020842
4287     if (isOldSorting) {
4288       BRepGProp::LinearProperties(S, GPr);
4289     }
4290     else {
4291       if (S.ShapeType() == TopAbs_EDGE || S.ShapeType() == TopAbs_WIRE) {
4292         BRepGProp::LinearProperties(S, GPr);
4293       }
4294       else if (S.ShapeType() == TopAbs_FACE || S.ShapeType() == TopAbs_SHELL) {
4295         BRepGProp::SurfaceProperties(S, GPr);
4296       }
4297       else {
4298         BRepGProp::VolumeProperties(S, GPr);
4299       }
4300     }
4301     // END: fix for Mantis issue 0020842
4302     GPoint = GPr.CentreOfMass();
4303     Len = GPr.Mass();
4304   }
4305
4306   double dMidXYZ = GPoint.X() * 999.0 + GPoint.Y() * 99.0 + GPoint.Z() * 0.9;
4307   return std::make_pair(dMidXYZ, Len);
4308 }
4309
4310 //=======================================================================
4311 //function : CompareShapes::operator()
4312 //purpose  : used by std::sort(), called from SortShapes()
4313 //=======================================================================
4314 bool GEOMImpl_IShapesOperations::CompareShapes::operator()(const TopoDS_Shape& theShape1,
4315                                                            const TopoDS_Shape& theShape2)
4316 {
4317   if (!myMap.IsBound(theShape1)) {
4318     myMap.Bind(theShape1, ShapeToDouble(theShape1, myIsOldSorting));
4319   }
4320
4321   if (!myMap.IsBound(theShape2)) {
4322     myMap.Bind(theShape2, ShapeToDouble(theShape2, myIsOldSorting));
4323   }
4324
4325   std::pair<double, double> val1 = myMap.Find(theShape1);
4326   std::pair<double, double> val2 = myMap.Find(theShape2);
4327
4328   double tol = Precision::Confusion();
4329   bool exchange = Standard_False;
4330
4331   double dMidXYZ = val1.first - val2.first;
4332   if (dMidXYZ >= tol) {
4333     exchange = Standard_True;
4334   }
4335   else if (Abs(dMidXYZ) < tol) {
4336     double dLength = val1.second - val2.second;
4337     if (dLength >= tol) {
4338       exchange = Standard_True;
4339     }
4340     else if (Abs(dLength) < tol && theShape1.ShapeType() <= TopAbs_FACE) {
4341       // PAL17233
4342       // equal values possible on shapes such as two halves of a sphere and
4343       // a membrane inside the sphere
4344       Bnd_Box box1,box2;
4345       BRepBndLib::Add(theShape1, box1);
4346       if (!box1.IsVoid()) {
4347         BRepBndLib::Add(theShape2, box2);
4348         Standard_Real dSquareExtent = box1.SquareExtent() - box2.SquareExtent();
4349         if (dSquareExtent >= tol) {
4350           exchange = Standard_True;
4351         }
4352         else if (Abs(dSquareExtent) < tol) {
4353           Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, val1, val2;
4354           box1.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
4355           val1 = (aXmin+aXmax)*999.0 + (aYmin+aYmax)*99.0 + (aZmin+aZmax)*0.9;
4356           box2.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
4357           val2 = (aXmin+aXmax)*999.0 + (aYmin+aYmax)*99.0 + (aZmin+aZmax)*0.9;
4358           if ((val1 - val2) >= tol) {
4359             exchange = Standard_True;
4360           }
4361         }
4362       }
4363     }
4364   }
4365
4366   //return val1 < val2;
4367   return !exchange;
4368 }
4369
4370 //=======================================================================
4371 //function : SortShapes
4372 //purpose  :
4373 //=======================================================================
4374 void GEOMImpl_IShapesOperations::SortShapes(TopTools_ListOfShape& SL,
4375                                             const Standard_Boolean isOldSorting)
4376 {
4377 #ifdef STD_SORT_ALGO
4378   std::vector<TopoDS_Shape> aShapesVec;
4379   aShapesVec.reserve(SL.Extent());
4380
4381   TopTools_ListIteratorOfListOfShape it (SL);
4382   for (; it.More(); it.Next()) {
4383     aShapesVec.push_back(it.Value());
4384   }
4385   SL.Clear();
4386
4387   CompareShapes shComp (isOldSorting);
4388   std::stable_sort(aShapesVec.begin(), aShapesVec.end(), shComp);
4389   //std::sort(aShapesVec.begin(), aShapesVec.end(), shComp);
4390
4391   std::vector<TopoDS_Shape>::const_iterator anIter = aShapesVec.begin();
4392   for (; anIter != aShapesVec.end(); ++anIter) {
4393     SL.Append(*anIter);
4394   }
4395 #else
4396   // old implementation
4397   Standard_Integer MaxShapes = SL.Extent();
4398   TopTools_Array1OfShape  aShapes (1,MaxShapes);
4399   TColStd_Array1OfInteger OrderInd(1,MaxShapes);
4400   TColStd_Array1OfReal    MidXYZ  (1,MaxShapes); //X,Y,Z;
4401   TColStd_Array1OfReal    Length  (1,MaxShapes); //X,Y,Z;
4402
4403   // Computing of CentreOfMass
4404   Standard_Integer Index;
4405   GProp_GProps GPr;
4406   gp_Pnt GPoint;
4407   TopTools_ListIteratorOfListOfShape it(SL);
4408   for (Index=1;  it.More();  Index++)
4409   {
4410     TopoDS_Shape S = it.Value();
4411     SL.Remove( it ); // == it.Next()
4412     aShapes(Index) = S;
4413     OrderInd.SetValue (Index, Index);
4414     if (S.ShapeType() == TopAbs_VERTEX) {
4415       GPoint = BRep_Tool::Pnt( TopoDS::Vertex( S ));
4416       Length.SetValue( Index, (Standard_Real) S.Orientation());
4417     }
4418     else {
4419       // BEGIN: fix for Mantis issue 0020842
4420       if (isOldSorting) {
4421         BRepGProp::LinearProperties (S, GPr);
4422       }
4423       else {
4424         if (S.ShapeType() == TopAbs_EDGE || S.ShapeType() == TopAbs_WIRE) {
4425           BRepGProp::LinearProperties (S, GPr);
4426         }
4427         else if (S.ShapeType() == TopAbs_FACE || S.ShapeType() == TopAbs_SHELL) {
4428           BRepGProp::SurfaceProperties(S, GPr);
4429         }
4430         else {
4431           BRepGProp::VolumeProperties(S, GPr);
4432         }
4433       }
4434       // END: fix for Mantis issue 0020842
4435       GPoint = GPr.CentreOfMass();
4436       Length.SetValue(Index, GPr.Mass());
4437     }
4438     MidXYZ.SetValue(Index, GPoint.X()*999.0 + GPoint.Y()*99.0 + GPoint.Z()*0.9);
4439     //cout << Index << " L: " << Length(Index) << "CG: " << MidXYZ(Index) << endl;
4440   }
4441
4442   // Sorting
4443   Standard_Integer aTemp;
4444   Standard_Boolean exchange, Sort = Standard_True;
4445   Standard_Real    tol = Precision::Confusion();
4446   while (Sort)
4447   {
4448     Sort = Standard_False;
4449     for (Index=1; Index < MaxShapes; Index++)
4450     {
4451       exchange = Standard_False;
4452       Standard_Real dMidXYZ = MidXYZ(OrderInd(Index)) - MidXYZ(OrderInd(Index+1));
4453       Standard_Real dLength = Length(OrderInd(Index)) - Length(OrderInd(Index+1));
4454       if ( dMidXYZ >= tol ) {
4455 //         cout << "MidXYZ: " << MidXYZ(OrderInd(Index))<< " > " <<MidXYZ(OrderInd(Index+1))
4456 //              << " d: " << dMidXYZ << endl;
4457         exchange = Standard_True;
4458       }
4459       else if ( Abs(dMidXYZ) < tol && dLength >= tol ) {
4460 //         cout << "Length: " << Length(OrderInd(Index))<< " > " <<Length(OrderInd(Index+1))
4461 //              << " d: " << dLength << endl;
4462         exchange = Standard_True;
4463       }
4464       else if ( Abs(dMidXYZ) < tol && Abs(dLength) < tol &&
4465                 aShapes(OrderInd(Index)).ShapeType() <= TopAbs_FACE) {
4466         // PAL17233
4467         // equal values possible on shapes such as two halves of a sphere and
4468         // a membrane inside the sphere
4469         Bnd_Box box1,box2;
4470         BRepBndLib::Add( aShapes( OrderInd(Index) ), box1 );
4471         if ( box1.IsVoid() ) continue;
4472         BRepBndLib::Add( aShapes( OrderInd(Index+1) ), box2 );
4473         Standard_Real dSquareExtent = box1.SquareExtent() - box2.SquareExtent();
4474         if ( dSquareExtent >= tol ) {
4475 //           cout << "SquareExtent: " << box1.SquareExtent()<<" > "<<box2.SquareExtent() << endl;
4476           exchange = Standard_True;
4477         }
4478         else if ( Abs(dSquareExtent) < tol ) {
4479           Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, val1, val2;
4480           box1.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
4481           val1 = (aXmin+aXmax)*999 + (aYmin+aYmax)*99 + (aZmin+aZmax)*0.9;
4482           box2.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
4483           val2 = (aXmin+aXmax)*999 + (aYmin+aYmax)*99 + (aZmin+aZmax)*0.9;
4484           //exchange = val1 > val2;
4485           if ((val1 - val2) >= tol) {
4486             exchange = Standard_True;
4487           }
4488           //cout << "box: " << val1<<" > "<<val2 << endl;
4489         }
4490       }
4491
4492       if (exchange)
4493       {
4494 //         cout << "exchange " << Index << " & " << Index+1 << endl;
4495         aTemp = OrderInd(Index);
4496         OrderInd(Index) = OrderInd(Index+1);
4497         OrderInd(Index+1) = aTemp;
4498         Sort = Standard_True;
4499       }
4500     }
4501   }
4502
4503   for (Index=1; Index <= MaxShapes; Index++)
4504     SL.Append( aShapes( OrderInd(Index) ));
4505 #endif
4506 }
4507
4508 //=======================================================================
4509 //function : CompsolidToCompound
4510 //purpose  :
4511 //=======================================================================
4512 TopoDS_Shape GEOMImpl_IShapesOperations::CompsolidToCompound (const TopoDS_Shape& theCompsolid)
4513 {
4514   if (theCompsolid.ShapeType() != TopAbs_COMPSOLID) {
4515     return theCompsolid;
4516   }
4517
4518   TopoDS_Compound aCompound;
4519   BRep_Builder B;
4520   B.MakeCompound(aCompound);
4521
4522   TopTools_MapOfShape mapShape;
4523   TopoDS_Iterator It (theCompsolid, Standard_True, Standard_True);
4524
4525   for (; It.More(); It.Next()) {
4526     TopoDS_Shape aShape_i = It.Value();
4527     if (mapShape.Add(aShape_i)) {
4528       B.Add(aCompound, aShape_i);
4529     }
4530   }
4531
4532   return aCompound;
4533 }
4534
4535 //=======================================================================
4536 //function : CheckTriangulation
4537 //purpose  :
4538 //=======================================================================
4539 bool GEOMImpl_IShapesOperations::CheckTriangulation (const TopoDS_Shape& aShape)
4540 {
4541   bool isTriangulation = true;
4542
4543   TopExp_Explorer exp (aShape, TopAbs_FACE);
4544   if (exp.More())
4545   {
4546     TopLoc_Location aTopLoc;
4547     Handle(Poly_Triangulation) aTRF;
4548     aTRF = BRep_Tool::Triangulation(TopoDS::Face(exp.Current()), aTopLoc);
4549     if (aTRF.IsNull()) {
4550       isTriangulation = false;
4551     }
4552   }
4553   else // no faces, try edges
4554   {
4555     TopExp_Explorer expe (aShape, TopAbs_EDGE);
4556     if (!expe.More()) {
4557       return false;
4558     }
4559     TopLoc_Location aLoc;
4560     Handle(Poly_Polygon3D) aPE = BRep_Tool::Polygon3D(TopoDS::Edge(expe.Current()), aLoc);
4561     if (aPE.IsNull()) {
4562       isTriangulation = false;
4563     }
4564   }
4565
4566   if (!isTriangulation) {
4567     // calculate deflection
4568     Standard_Real aDeviationCoefficient = 0.001;
4569
4570     Bnd_Box B;
4571     BRepBndLib::Add(aShape, B);
4572     Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
4573     B.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
4574
4575     Standard_Real dx = aXmax - aXmin, dy = aYmax - aYmin, dz = aZmax - aZmin;
4576     Standard_Real aDeflection = Max(Max(dx, dy), dz) * aDeviationCoefficient * 4;
4577     Standard_Real aHLRAngle = 0.349066;
4578
4579     BRepMesh_IncrementalMesh Inc (aShape, aDeflection, Standard_False, aHLRAngle);
4580   }
4581
4582   return true;
4583 }
4584
4585 #define MAX_TOLERANCE 1.e-7
4586
4587 //=======================================================================
4588 //function : isSameEdge
4589 //purpose  : Returns True if two edges coincide
4590 //=======================================================================
4591 static bool isSameEdge(const TopoDS_Edge& theEdge1, const TopoDS_Edge& theEdge2)
4592 {
4593   TopoDS_Vertex V11, V12, V21, V22;
4594   TopExp::Vertices(theEdge1, V11, V12);
4595   TopExp::Vertices(theEdge2, V21, V22);
4596   gp_Pnt P11 = BRep_Tool::Pnt(V11);
4597   gp_Pnt P12 = BRep_Tool::Pnt(V12);
4598   gp_Pnt P21 = BRep_Tool::Pnt(V21);
4599   gp_Pnt P22 = BRep_Tool::Pnt(V22);
4600   bool coincide = false;
4601
4602   //Check that ends of edges coincide
4603   if(P11.Distance(P21) <= MAX_TOLERANCE) {
4604     if(P12.Distance(P22) <= MAX_TOLERANCE) coincide =  true;
4605   }
4606   else if(P11.Distance(P22) <= MAX_TOLERANCE) {
4607     if(P12.Distance(P21) <= MAX_TOLERANCE) coincide = true;
4608   }
4609
4610   if(!coincide) return false;
4611
4612   if (BRep_Tool::Degenerated(theEdge1))
4613     if (BRep_Tool::Degenerated(theEdge2)) return true;
4614     else return false;
4615   else
4616     if (BRep_Tool::Degenerated(theEdge2)) return false;
4617
4618   double U11, U12, U21, U22;
4619   Handle(Geom_Curve) C1 = BRep_Tool::Curve(theEdge1, U11, U12);
4620   Handle(Geom_Curve) C2 = BRep_Tool::Curve(theEdge2, U21, U22);
4621   if(C1->DynamicType() == C2->DynamicType()) return true;
4622
4623   //Check that both edges has the same geometry
4624   double range = U12-U11;
4625   double U = U11+ range/3.0;
4626   gp_Pnt P1 = C1->Value(U);     //Compute a point on one third of the edge's length
4627   U = U11+range*2.0/3.0;
4628   gp_Pnt P2 = C1->Value(U);     //Compute a point on two thirds of the edge's length
4629
4630   if(!GeomLib_Tool::Parameter(C2, P1, MAX_TOLERANCE, U) ||  U < U21 || U > U22)
4631     return false;
4632
4633   if(P1.Distance(C2->Value(U)) > MAX_TOLERANCE) return false;
4634
4635   if(!GeomLib_Tool::Parameter(C2, P2, MAX_TOLERANCE, U) || U < U21 || U > U22)
4636     return false;
4637
4638   if(P2.Distance(C2->Value(U)) > MAX_TOLERANCE) return false;
4639
4640   return true;
4641 }
4642
4643 #include <TopoDS_TShape.hxx>
4644 //=======================================================================
4645 //function : isSameFace
4646 //purpose  : Returns True if two faces coincide
4647 //=======================================================================
4648 static bool isSameFace(const TopoDS_Face& theFace1, const TopoDS_Face& theFace2)
4649 {
4650   TopExp_Explorer E(theFace1, TopAbs_EDGE);
4651   TopTools_ListOfShape LS1, LS2;
4652   for(; E.More(); E.Next()) LS1.Append(E.Current());
4653
4654   E.Init(theFace2, TopAbs_EDGE);
4655   for(; E.More(); E.Next()) LS2.Append(E.Current());
4656
4657   //Compare the number of edges in the faces
4658   if(LS1.Extent() != LS2.Extent()) return false;
4659
4660   double aMin = RealFirst(), aMax = RealLast();
4661   double xminB1=aMax, yminB1=aMax, zminB1=aMax, xminB2=aMax, yminB2=aMax, zminB2=aMax;
4662   double xmaxB1=aMin, ymaxB1=aMin, zmaxB1=aMin, xmaxB2=aMin, ymaxB2=aMin, zmaxB2=aMin;
4663
4664   for(E.Init(theFace1, TopAbs_VERTEX); E.More(); E.Next()) {
4665     gp_Pnt P = BRep_Tool::Pnt(TopoDS::Vertex(E.Current()));
4666     if(P.X() < xminB1) xminB1 = P.X();
4667     if(P.Y() < yminB1) yminB1 = P.Y();
4668     if(P.Z() < zminB1) zminB1 = P.Z();
4669     if(P.X() > xmaxB1) xmaxB1 = P.X();
4670     if(P.Y() > ymaxB1) ymaxB1 = P.Y();
4671     if(P.Z() > zmaxB1) zmaxB1 = P.Z();
4672   }
4673
4674   for(E.Init(theFace2, TopAbs_VERTEX); E.More(); E.Next()) {
4675     gp_Pnt P = BRep_Tool::Pnt(TopoDS::Vertex(E.Current()));
4676     if(P.X() < xminB2) xminB2 = P.X();
4677     if(P.Y() < yminB2) yminB2 = P.Y();
4678     if(P.Z() < zminB2) zminB2 = P.Z();
4679     if(P.X() > xmaxB2) xmaxB2 = P.X();
4680     if(P.Y() > ymaxB2) ymaxB2 = P.Y();
4681     if(P.Z() > zmaxB2) zmaxB2 = P.Z();
4682   }
4683
4684   //Compare the bounding boxes of both faces
4685   if(gp_Pnt(xminB1, yminB1, zminB1).Distance(gp_Pnt(xminB2, yminB2, zminB2)) > MAX_TOLERANCE)
4686     return false;
4687
4688   if(gp_Pnt(xmaxB1, ymaxB1, zmaxB1).Distance(gp_Pnt(xmaxB2, ymaxB2, zmaxB2)) > MAX_TOLERANCE)
4689     return false;
4690
4691   Handle(Geom_Surface) S1 = BRep_Tool::Surface(theFace1);
4692   Handle(Geom_Surface) S2 = BRep_Tool::Surface(theFace2);
4693
4694   //Check if there a coincidence of two surfaces at least in two points
4695   double U11, U12, V11, V12, U21, U22, V21, V22;
4696   BRepTools::UVBounds(theFace1, U11, U12, V11, V12);
4697   BRepTools::UVBounds(theFace2, U21, U22, V21, V22);
4698
4699   double rangeU = U12-U11;
4700   double rangeV = V12-V11;
4701   double U = U11 + rangeU/3.0;
4702   double V = V11 + rangeV/3.0;
4703   gp_Pnt P1 = S1->Value(U, V);
4704   U = U11+rangeU*2.0/3.0;
4705   V = V11+rangeV*2.0/3.0;
4706   gp_Pnt P2 = S1->Value(U, V);
4707
4708   if (!GeomLib_Tool::Parameters(S2, P1, MAX_TOLERANCE, U, V) || U < U21 || U > U22 || V < V21 || V > V22)
4709     return false;
4710
4711   if (P1.Distance(S2->Value(U,V)) > MAX_TOLERANCE) return false;
4712
4713   if (!GeomLib_Tool::Parameters(S2, P2, MAX_TOLERANCE, U, V) || U < U21 || U > U22 || V < V21 || V > V22)
4714     return false;
4715
4716   if (P2.Distance(S2->Value(U, V)) > MAX_TOLERANCE) return false;
4717
4718   //Check that each edge of the Face1 has a counterpart in the Face2
4719   TopTools_MapOfOrientedShape aMap;
4720   TopTools_ListIteratorOfListOfShape LSI1(LS1);
4721   for(; LSI1.More(); LSI1.Next()) {
4722     TopoDS_Edge E = TopoDS::Edge(LSI1.Value());
4723     bool isFound = false;
4724     TopTools_ListIteratorOfListOfShape LSI2(LS2);
4725     for(; LSI2.More(); LSI2.Next()) {
4726       TopoDS_Shape aValue = LSI2.Value();
4727       if(aMap.Contains(aValue)) continue; //To avoid checking already found edge several times
4728       if(isSameEdge(E, TopoDS::Edge(aValue))) {
4729         aMap.Add(aValue);
4730         isFound = true;
4731         break;
4732       }
4733     }
4734     if(!isFound) return false;
4735   }
4736
4737   return true;
4738 }
4739
4740 //=======================================================================
4741 //function : isSameSolid
4742 //purpose  : Returns True if two solids coincide
4743 //=======================================================================
4744 bool isSameSolid(const TopoDS_Solid& theSolid1, const TopoDS_Solid& theSolid2)
4745 {
4746   TopExp_Explorer E(theSolid1, TopAbs_FACE);
4747   TopTools_ListOfShape LS1, LS2;
4748   for(; E.More(); E.Next()) LS1.Append(E.Current());
4749   E.Init(theSolid2, TopAbs_FACE);
4750   for(; E.More(); E.Next()) LS2.Append(E.Current());
4751
4752   if(LS1.Extent() != LS2.Extent()) return false;
4753
4754   double aMin = RealFirst(), aMax = RealLast();
4755   double xminB1=aMax, yminB1=aMax, zminB1=aMax, xminB2=aMax, yminB2=aMax, zminB2=aMax;
4756   double xmaxB1=aMin, ymaxB1=aMin, zmaxB1=aMin, xmaxB2=aMin, ymaxB2=aMin, zmaxB2=aMin;
4757
4758   for(E.Init(theSolid1, TopAbs_VERTEX); E.More(); E.Next()) {
4759     gp_Pnt P = BRep_Tool::Pnt(TopoDS::Vertex(E.Current()));
4760     if(P.X() < xminB1) xminB1 = P.X();
4761     if(P.Y() < yminB1) yminB1 = P.Y();
4762     if(P.Z() < zminB1) zminB1 = P.Z();
4763     if(P.X() > xmaxB1) xmaxB1 = P.X();
4764     if(P.Y() > ymaxB1) ymaxB1 = P.Y();
4765     if(P.Z() > zmaxB1) zmaxB1 = P.Z();
4766   }
4767
4768   for(E.Init(theSolid2, TopAbs_VERTEX); E.More(); E.Next()) {
4769     gp_Pnt P = BRep_Tool::Pnt(TopoDS::Vertex(E.Current()));
4770     if(P.X() < xminB2) xminB2 = P.X();
4771     if(P.Y() < yminB2) yminB2 = P.Y();
4772     if(P.Z() < zminB2) zminB2 = P.Z();
4773     if(P.X() > xmaxB2) xmaxB2 = P.X();
4774     if(P.Y() > ymaxB2) ymaxB2 = P.Y();
4775     if(P.Z() > zmaxB2) zmaxB2 = P.Z();
4776   }
4777
4778   //Compare the bounding boxes of both solids
4779   if(gp_Pnt(xminB1, yminB1, zminB1).Distance(gp_Pnt(xminB2, yminB2, zminB2)) > MAX_TOLERANCE)
4780     return false;
4781
4782   if(gp_Pnt(xmaxB1, ymaxB1, zmaxB1).Distance(gp_Pnt(xmaxB2, ymaxB2, zmaxB2)) > MAX_TOLERANCE)
4783     return false;
4784
4785   //Check that each face of the Solid1 has a counterpart in the Solid2
4786   TopTools_MapOfOrientedShape aMap;
4787   TopTools_ListIteratorOfListOfShape LSI1(LS1);
4788   for(; LSI1.More(); LSI1.Next()) {
4789     TopoDS_Face F = TopoDS::Face(LSI1.Value());
4790     bool isFound = false;
4791     TopTools_ListIteratorOfListOfShape LSI2(LS2);
4792     for(; LSI2.More(); LSI2.Next()) {
4793       if(aMap.Contains(LSI2.Value())) continue; //To avoid checking already found faces several times
4794       if(isSameFace(F, TopoDS::Face(LSI2.Value()))) {
4795         aMap.Add(LSI2.Value());
4796         isFound = true;
4797         break;
4798       }
4799     }
4800     if(!isFound) return false;
4801   }
4802
4803   return true;
4804 }
4805
4806 //=======================================================================
4807 //function : GetSame
4808 //purpose  :
4809 //=======================================================================
4810 Handle(GEOM_Object) GEOMImpl_IShapesOperations::GetSame(const Handle(GEOM_Object)& theShapeWhere,
4811                                                         const Handle(GEOM_Object)& theShapeWhat)
4812 {
4813   SetErrorCode(KO);
4814   if (theShapeWhere.IsNull() || theShapeWhat.IsNull()) return NULL;
4815
4816   TopoDS_Shape aWhere = theShapeWhere->GetValue();
4817   TopoDS_Shape aWhat  = theShapeWhat->GetValue();
4818
4819   if (aWhere.IsNull() || aWhat.IsNull()) return NULL;
4820
4821   int anIndex = -1;
4822   bool isFound = false;
4823   TopoDS_Shape aSubShape;
4824   TopTools_MapOfShape aMap;
4825
4826   if (aWhat.ShapeType() == TopAbs_COMPOUND || aWhat.ShapeType() == TopAbs_COMPSOLID) {
4827     TopoDS_Iterator It (aWhat, Standard_True, Standard_True);
4828     if (It.More()) aWhat = It.Value();
4829     It.Next();
4830     if (It.More()) {
4831       SetErrorCode("Compounds of two or more shapes are not allowed for aWhat argument");
4832       return NULL;
4833     }
4834   }
4835
4836   switch (aWhat.ShapeType()) {
4837     case TopAbs_VERTEX: {
4838       gp_Pnt P = BRep_Tool::Pnt(TopoDS::Vertex(aWhat));
4839       TopExp_Explorer E(aWhere, TopAbs_VERTEX);
4840       for(; E.More(); E.Next()) {
4841         if(!aMap.Add(E.Current())) continue;
4842         gp_Pnt P2 = BRep_Tool::Pnt(TopoDS::Vertex(E.Current()));
4843         if(P.Distance(P2) <= MAX_TOLERANCE) {
4844           isFound = true;
4845           aSubShape = E.Current();
4846           break;
4847         }
4848       }
4849       break;
4850                         }
4851     case TopAbs_EDGE: {
4852       TopoDS_Edge anEdge = TopoDS::Edge(aWhat);
4853       TopExp_Explorer E(aWhere, TopAbs_EDGE);
4854       for(; E.More(); E.Next()) {
4855         if(!aMap.Add(E.Current())) continue;
4856         if(isSameEdge(anEdge, TopoDS::Edge(E.Current()))) {
4857           aSubShape = E.Current();
4858           isFound = true;
4859           break;
4860         }
4861       }
4862       break;
4863                       }
4864     case TopAbs_FACE: {
4865       TopoDS_Face aFace = TopoDS::Face(aWhat);
4866       TopExp_Explorer E(aWhere, TopAbs_FACE);
4867       for(; E.More(); E.Next()) {
4868         if(!aMap.Add(E.Current())) continue;
4869         if(isSameFace(aFace, TopoDS::Face(E.Current()))) {
4870           aSubShape = E.Current();
4871           isFound = true;
4872           break;
4873         }
4874       }
4875       break;
4876                       }
4877     case TopAbs_SOLID: {
4878       TopoDS_Solid aSolid = TopoDS::Solid(aWhat);
4879       TopExp_Explorer E(aWhere, TopAbs_SOLID);
4880       for(; E.More(); E.Next()) {
4881         if(!aMap.Add(E.Current())) continue;
4882         if(isSameSolid(aSolid, TopoDS::Solid(E.Current()))) {
4883           aSubShape = E.Current();
4884           isFound = true;
4885           break;
4886         }
4887       }
4888       break;
4889                        }
4890     default:
4891       return NULL;
4892   }
4893
4894   if (isFound) {
4895     TopTools_IndexedMapOfShape anIndices;
4896     TopExp::MapShapes(aWhere, anIndices);
4897     if (anIndices.Contains(aSubShape))
4898       anIndex = anIndices.FindIndex(aSubShape);
4899   }
4900
4901   if (anIndex < 0) return NULL;
4902
4903   Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
4904
4905   anArray->SetValue(1, anIndex);
4906
4907   Handle(GEOM_Object) aResult = GetEngine()->AddSubShape(theShapeWhere, anArray);
4908   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
4909
4910   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetSame("
4911     << theShapeWhere << ", " << theShapeWhat << ")";
4912
4913   SetErrorCode(OK);
4914
4915   return aResult;
4916 }
4917
4918
4919 //=======================================================================
4920 //function : GetSameIDs
4921 //purpose  :
4922 //=======================================================================
4923 Handle(TColStd_HSequenceOfInteger) GEOMImpl_IShapesOperations::GetSameIDs(const Handle(GEOM_Object)& theShapeWhere,
4924                                                         const Handle(GEOM_Object)& theShapeWhat)
4925 {
4926   SetErrorCode(KO);
4927   if (theShapeWhere.IsNull() || theShapeWhat.IsNull()) return NULL;
4928
4929   TopoDS_Shape aWhere = theShapeWhere->GetValue();
4930   TopoDS_Shape aWhat  = theShapeWhat->GetValue();
4931
4932   if (aWhere.IsNull() || aWhat.IsNull()) return NULL;
4933
4934   int anIndex = -1;
4935   bool isFound = false;
4936   TopTools_ListOfShape listShape;
4937   TopTools_MapOfShape aMap;
4938
4939   if (aWhat.ShapeType() == TopAbs_COMPOUND || aWhat.ShapeType() == TopAbs_COMPSOLID) {
4940     TopoDS_Iterator It (aWhat, Standard_True, Standard_True);
4941     if (It.More()) aWhat = It.Value();
4942     It.Next();
4943     if (It.More()) {
4944       SetErrorCode("Compounds of two or more shapes are not allowed for aWhat argument");
4945       return NULL;
4946     }
4947   }
4948
4949   switch (aWhat.ShapeType()) {
4950     case TopAbs_VERTEX: {
4951       gp_Pnt P = BRep_Tool::Pnt(TopoDS::Vertex(aWhat));
4952       TopExp_Explorer E(aWhere, TopAbs_VERTEX);
4953       for(; E.More(); E.Next()) {
4954         if(!aMap.Add(E.Current())) continue;
4955         gp_Pnt P2 = BRep_Tool::Pnt(TopoDS::Vertex(E.Current()));
4956         if(P.Distance(P2) <= MAX_TOLERANCE) {
4957           listShape.Append(E.Current());
4958         }
4959       }
4960       break;
4961                         }
4962     case TopAbs_EDGE: {
4963       TopoDS_Edge anEdge = TopoDS::Edge(aWhat);
4964       TopExp_Explorer E(aWhere, TopAbs_EDGE);
4965       for(; E.More(); E.Next()) {
4966         if(!aMap.Add(E.Current())) continue;
4967         if(isSameEdge(anEdge, TopoDS::Edge(E.Current()))) {
4968           listShape.Append(E.Current());
4969         }
4970       }
4971       break;
4972                       }
4973     case TopAbs_FACE: {
4974       TopoDS_Face aFace = TopoDS::Face(aWhat);
4975       TopExp_Explorer E(aWhere, TopAbs_FACE);
4976       for(; E.More(); E.Next()) {
4977         if(!aMap.Add(E.Current())) continue;
4978         if(isSameFace(aFace, TopoDS::Face(E.Current()))) {
4979           listShape.Append(E.Current());
4980         }
4981       }
4982       break;
4983                       }
4984     case TopAbs_SOLID: {
4985       TopoDS_Solid aSolid = TopoDS::Solid(aWhat);
4986       TopExp_Explorer E(aWhere, TopAbs_SOLID);
4987       for(; E.More(); E.Next()) {
4988         if(!aMap.Add(E.Current())) continue;
4989         if(isSameSolid(aSolid, TopoDS::Solid(E.Current()))) {
4990           listShape.Append(E.Current());
4991         }
4992       }
4993       break;
4994                        }
4995     default:
4996       return NULL;
4997   }
4998
4999   if ( !listShape.IsEmpty() ) {
5000     TopTools_IndexedMapOfShape anIndices;
5001     TopExp::MapShapes(aWhere, anIndices);
5002     TopTools_ListIteratorOfListOfShape itSub (listShape);
5003     Handle(TColStd_HSequenceOfInteger) aSeq = new TColStd_HSequenceOfInteger;
5004     for (; itSub.More(); itSub.Next()) {
5005       if (anIndices.Contains(itSub.Value()))
5006         aSeq->Append(anIndices.FindIndex(itSub.Value()));
5007     }
5008     SetErrorCode(OK);
5009     // The GetSameIDs() doesn't change object so no new function is required.
5010     Handle(GEOM_Function) aFunction = GEOM::GetCreatedLast(theShapeWhere,theShapeWhat)->GetLastFunction();
5011
5012   // Make a Python command
5013   GEOM::TPythonDump(aFunction)
5014     << "listSameIDs = geompy.GetSameIDs("
5015     << theShapeWhere << ", "
5016     << theShapeWhat << ")";
5017     return aSeq;
5018   } else {
5019     SetErrorCode(NOT_FOUND_ANY);
5020     return NULL; 
5021   }
5022 }