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