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