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