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