bc1a0ba73fc6fc0103c0555efedf2aa678cc335c
[occt.git] / src / BOPTools / BOPTools_AlgoTools.cxx
1 // Created by: Peter KURNEV
2 // Copyright (c) 2010-2012 OPEN CASCADE SAS
3 // Copyright (c) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE
4 // Copyright (c) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, CEDRAT,
5 //                         EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 //
7 // The content of this file is subject to the Open CASCADE Technology Public
8 // License Version 6.5 (the "License"). You may not use the content of this file
9 // except in compliance with the License. Please obtain a copy of the License
10 // at http://www.opencascade.org and read it completely before using this file.
11 //
12 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
13 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
14 //
15 // The Original Code and all software distributed under the License is
16 // distributed on an "AS IS" basis, without warranty of any kind, and the
17 // Initial Developer hereby disclaims all such warranties, including without
18 // limitation, any warranties of merchantability, fitness for a particular
19 // purpose or non-infringement. Please see the License for the specific terms
20 // and conditions governing the rights and limitations under the License.
21
22
23 #include <BOPTools_AlgoTools.ixx>
24 //
25 #include <Precision.hxx>
26 //
27 #include <gp_Pnt.hxx>
28 #include <gp_XYZ.hxx>
29 #include <gp_Pnt2d.hxx>
30 #include <gp_Cylinder.hxx>
31 #include <gp_Cone.hxx>
32 #include <gp_Sphere.hxx>
33 #include <gp_Torus.hxx>
34 //
35 #include <Geom2d_Curve.hxx>
36 #include <Geom_Surface.hxx>
37 #include <Geom_Plane.hxx>
38 #include <Geom_TrimmedCurve.hxx>
39 #include <Geom_Curve.hxx>
40 #include <GeomAPI_ProjectPointOnSurf.hxx>
41 #include <Geom2dInt_Geom2dCurveTool.hxx>
42 #include <GeomAdaptor_Surface.hxx>
43 //
44 #include <TopAbs_Orientation.hxx>
45 //
46 #include <TopoDS_Compound.hxx>
47 #include <TopoDS_CompSolid.hxx>
48 #include <TopoDS_Solid.hxx>
49 #include <TopoDS_Shell.hxx>
50 #include <TopoDS_Wire.hxx>
51 //
52 #include <BRep_Builder.hxx>
53 #include <BRep_Tool.hxx>
54 #include <BRepLib.hxx>
55 #include <BRepAdaptor_Curve2d.hxx>
56 #include <BRepClass3d_SolidClassifier.hxx>
57 #include <TopExp.hxx>
58 #include <TopExp_Explorer.hxx>
59 //
60 #include <IntTools_Tools.hxx>
61 //
62 #include <BOPTools.hxx>
63 #include <BOPTools_CoupleOfShape.hxx>
64 #include <BOPTools_ListOfCoupleOfShape.hxx>
65 #include <BOPTools_AlgoTools2D.hxx>
66 #include <BOPTools_AlgoTools3D.hxx>
67 //
68 #include <BOPCol_IndexedMapOfShape.hxx>
69 #include <BOPCol_MapOfShape.hxx>
70 //
71 #include <BOPInt_ShrunkRange.hxx>
72 //
73
74 static
75   Standard_Real AngleWithRef(const gp_Dir& theD1,
76                              const gp_Dir& theD2,
77                              const gp_Dir& theDRef);
78 static 
79   inline Standard_Real fsqrt(Standard_Real val);
80
81 static
82   Standard_Boolean FindFacePairs (const TopoDS_Edge& theE,
83                                   const BOPCol_ListOfShape& thLF,
84                                   BOPTools_ListOfCoupleOfShape& theLCFF,
85                                   Handle(BOPInt_Context)& theContext);
86 static
87   TopAbs_Orientation Orientation(const TopoDS_Edge& anE,
88                                  const TopoDS_Face& aF);
89
90 static 
91   void GetFaceDir(const TopoDS_Edge& aE,
92                   const TopoDS_Face& aF,
93                   const gp_Pnt& aP,
94                   const Standard_Real aT,
95                   const gp_Dir& aDTgt,
96                   gp_Dir& aDN,
97                   gp_Dir& aDB,
98                   Handle(BOPInt_Context)& theContext);
99 static 
100   Standard_Boolean FindPointInFace(const TopoDS_Edge& aE,
101                                    const TopoDS_Face& aF,
102                                    const gp_Pnt& aP,
103                                    const Standard_Real aT,
104                                    gp_Dir& aDB,
105                                    gp_Pnt& aPOut,
106                                    Handle(BOPInt_Context)& theContext,
107                                    const GeomAdaptor_Surface& aGAS);
108
109 //=======================================================================
110 // function: MakeConnexityBlocks
111 // purpose: 
112 //=======================================================================
113 void BOPTools_AlgoTools::MakeConnexityBlocks (const TopoDS_Shape& theS,
114                                               const TopAbs_ShapeEnum theType1,
115                                               const TopAbs_ShapeEnum theType2,
116                                               BOPCol_ListOfShape& theLCB)
117 {
118   Standard_Integer  aNbF, aNbAdd, aNbAdd1, i;
119   BRep_Builder aBB;
120   TopoDS_Compound aC;
121   TopoDS_Iterator aIt;
122   TopExp_Explorer aExp;
123   BOPCol_MapOfShape aMP;
124   BOPCol_IndexedMapOfShape aMCB, aMAdd, aMAdd1;
125   BOPCol_IndexedDataMapOfShapeListOfShape aMEF;
126   BOPCol_ListIteratorOfListOfShape aItLF;
127   //
128   // 1. aMEF
129   BOPTools::MapShapesAndAncestors(theS, theType1, theType2, aMEF);
130   //
131   // 2. aMCB
132   aIt.Initialize(theS);
133   for (; aIt.More(); aIt.Next()) {
134     const TopoDS_Shape& aF1=aIt.Value();
135     if (aMP.Contains(aF1)) {
136       continue;
137     }
138     //
139     aMCB.Clear();
140     aMAdd.Clear();
141     aMAdd.Add(aF1);
142     //
143     for(;;) {
144       aMAdd1.Clear();
145       //
146       aNbAdd = aMAdd.Extent();
147       for (i=1; i<=aNbAdd; ++i) {
148         const TopoDS_Shape& aF=aMAdd(i);
149         //
150         aExp.Init(aF, theType1);
151         for (; aExp.More(); aExp.Next()) {
152           const TopoDS_Shape& aE=aExp.Current();
153           //
154           const BOPCol_ListOfShape& aLF=aMEF.FindFromKey(aE);
155           aItLF.Initialize(aLF);
156           for (; aItLF.More(); aItLF.Next()) {
157             const TopoDS_Shape& aFx=aItLF.Value();
158             if (aFx.IsSame(aF)) {
159               continue;
160             }
161             if (aMCB.Contains(aFx)) {
162               continue;
163             }
164             aMAdd1.Add(aFx);
165           }
166         }//for (; aExp.More(); aExp.Next()){
167         aMCB.Add(aF);
168       }// for (i=1; i<=aNbAdd; ++i) {
169       //
170       aNbAdd1=aMAdd1.Extent();
171       if (!aNbAdd1) {
172         break;// ->make new CB from aMCB
173       }
174       //
175       aMAdd.Clear();
176       for (i=1; i<=aNbAdd1; ++i) {
177         const TopoDS_Shape& aFAdd = aMAdd1(i);
178         aMAdd.Add(aFAdd);
179       }
180     }//while(1) {
181     //
182     aNbF=aMCB.Extent();
183     if (aNbF) {
184       aBB.MakeCompound(aC);
185       //
186       for (i=1; i<=aNbF; ++i) {
187         const TopoDS_Shape& aF=aMCB(i);
188         aBB.Add(aC, aF);  
189         aMP.Add(aF);
190       }
191       theLCB.Append(aC);
192     }
193   }// for (; aIt.More(); aIt.Next()) 
194 }
195 //=======================================================================
196 // function: OrientFacesOnShell
197 // purpose: 
198 //=======================================================================
199 void BOPTools_AlgoTools::OrientFacesOnShell (TopoDS_Shape& aShell)
200 {
201   Standard_Boolean bIsProcessed1, bIsProcessed2;
202   Standard_Integer i, aNbE, aNbF, j;
203   TopAbs_Orientation anOrE1, anOrE2;
204   TopoDS_Face aF1x, aF2x;
205   TopoDS_Shape aShellNew;
206   BOPCol_IndexedDataMapOfShapeListOfShape aEFMap;
207   BOPCol_IndexedMapOfShape aProcessedFaces;
208   BRep_Builder aBB;
209   //
210   BOPTools_AlgoTools::MakeContainer(TopAbs_SHELL, aShellNew);
211   //
212   BOPTools::MapShapesAndAncestors(aShell, TopAbs_EDGE, TopAbs_FACE, aEFMap);
213   aNbE=aEFMap.Extent();
214   // 
215   // One seam edge  in aEFMap contains  2 equivalent faces.
216   for (i=1; i<=aNbE; ++i) {
217     BOPCol_ListOfShape& aLF=aEFMap.ChangeFromIndex(i);
218     aNbF=aLF.Extent();
219     if (aNbF>1) {
220       BOPCol_ListOfShape aLFTmp;
221       BOPCol_IndexedMapOfShape aFM;
222       //
223       BOPCol_ListIteratorOfListOfShape anIt(aLF);
224       for (; anIt.More(); anIt.Next()) {
225         const TopoDS_Shape& aF=anIt.Value();
226         if (!aFM.Contains(aF)) {
227           aFM.Add(aF);
228           aLFTmp.Append(aF);
229         }
230       }
231       aLF.Clear();
232       aLF=aLFTmp;
233     }
234   }
235   //
236   // Do
237   for (i=1; i<=aNbE; ++i) {
238     const TopoDS_Edge& aE=(*(TopoDS_Edge*)(&aEFMap.FindKey(i)));
239     if (BRep_Tool::Degenerated(aE)) {
240       continue;
241     }
242     //
243     const BOPCol_ListOfShape& aLF=aEFMap.FindFromIndex(i);
244     aNbF=aLF.Extent();
245     if (aNbF!=2) {
246       continue;
247     }
248     //
249     TopoDS_Face& aF1=(*(TopoDS_Face*)(&aLF.First()));
250     TopoDS_Face& aF2=(*(TopoDS_Face*)(&aLF.Last()));
251     //    
252     bIsProcessed1=aProcessedFaces.Contains(aF1);
253     bIsProcessed2=aProcessedFaces.Contains(aF2);
254     if (bIsProcessed1 && bIsProcessed2) {
255       continue;
256     }
257     
258     if (!bIsProcessed1 && !bIsProcessed2) {
259       aProcessedFaces.Add(aF1);
260       aBB.Add(aShellNew, aF1);
261       bIsProcessed1=!bIsProcessed1;
262     }
263     //
264     aF1x=aF1;
265     if (bIsProcessed1) {
266       j=aProcessedFaces.FindIndex(aF1);
267       aF1x=(*(TopoDS_Face*)(&aProcessedFaces.FindKey(j)));
268     }
269     //
270     aF2x=aF2;
271     if (bIsProcessed2) {
272       j=aProcessedFaces.FindIndex(aF2);
273       aF2x=(*(TopoDS_Face*)(&aProcessedFaces.FindKey(j)));
274     }
275     //
276     anOrE1=Orientation(aE, aF1x); 
277     anOrE2=Orientation(aE, aF2x);
278     //
279     if (bIsProcessed1 && !bIsProcessed2) {
280       if (anOrE1==anOrE2) {
281         if (!BRep_Tool::IsClosed(aE, aF1) &&
282             !BRep_Tool::IsClosed(aE, aF2)) {
283           aF2.Reverse();
284         }
285       }
286       aProcessedFaces.Add(aF2);
287       aBB.Add(aShellNew, aF2);
288     }
289     else if (!bIsProcessed1 && bIsProcessed2) {
290       if (anOrE1==anOrE2) {
291         if (!BRep_Tool::IsClosed(aE, aF1) &&
292             !BRep_Tool::IsClosed(aE, aF2)) {
293           aF1.Reverse();
294         }
295       }
296       aProcessedFaces.Add(aF1);
297       aBB.Add(aShellNew, aF1);
298     }
299   }
300   //
301   //
302   for (i=1; i<=aNbE; ++i) {
303     const TopoDS_Edge& aE=(*(TopoDS_Edge*)(&aEFMap.FindKey(i)));
304     if (BRep_Tool::Degenerated(aE)) {
305       continue;
306     }
307     //
308     const BOPCol_ListOfShape& aLF=aEFMap.FindFromIndex(i);
309     aNbF=aLF.Extent();
310     if (aNbF!=2) {
311       BOPCol_ListIteratorOfListOfShape anIt(aLF);
312       for(; anIt.More(); anIt.Next()) {
313         const TopoDS_Face& aF=(*(TopoDS_Face*)(&anIt.Value()));
314         if (!aProcessedFaces.Contains(aF)) {
315           aProcessedFaces.Add(aF);
316           aBB.Add(aShellNew, aF);
317         }
318       }
319     }
320   }
321   aShell=aShellNew;
322 }
323 //=======================================================================
324 //function : Orientation
325 //purpose  :
326 //=======================================================================
327 TopAbs_Orientation Orientation(const TopoDS_Edge& anE,
328                                const TopoDS_Face& aF)
329 {
330   TopAbs_Orientation anOr=TopAbs_INTERNAL;
331
332   TopExp_Explorer anExp;
333   anExp.Init(aF, TopAbs_EDGE);
334   for (; anExp.More(); anExp.Next()) {
335     const TopoDS_Edge& anEF1=(*(TopoDS_Edge*)(&anExp.Current()));
336     if (anEF1.IsSame(anE)) {
337       anOr=anEF1.Orientation();
338       break;
339     }
340   }
341   return anOr;
342 }
343
344 ////////////
345
346
347 //=======================================================================
348 // function: MakeConnexityBlock.
349 // purpose: 
350 //=======================================================================
351 void BOPTools_AlgoTools::MakeConnexityBlock (BOPCol_ListOfShape& theLFIn,
352                                              BOPCol_IndexedMapOfShape& theMEAvoid,
353                                              BOPCol_ListOfShape& theLCB,
354                                              const Handle(NCollection_BaseAllocator)& theAllocator)
355 {
356   Standard_Integer  aNbF, aNbAdd1, aNbAdd, i;
357   TopExp_Explorer aExp;
358   BOPCol_ListIteratorOfListOfShape aIt;
359   //
360   BOPCol_IndexedMapOfShape aMCB(100, theAllocator);
361   BOPCol_IndexedMapOfShape aMAdd(100, theAllocator);
362   BOPCol_IndexedMapOfShape aMAdd1(100, theAllocator);
363   BOPCol_IndexedDataMapOfShapeListOfShape aMEF(100, theAllocator);
364   //
365   // 1. aMEF
366   aNbF=theLFIn.Extent();
367   aIt.Initialize(theLFIn);
368   for (; aIt.More(); aIt.Next()) {
369     const TopoDS_Shape& aF=aIt.Value();      
370     BOPTools::MapShapesAndAncestors(aF, TopAbs_EDGE, TopAbs_FACE, aMEF);
371   }
372   //
373   // 2. aMCB
374   const TopoDS_Shape& aF1=theLFIn.First();
375   aMAdd.Add(aF1);
376   //
377   for(;;) {
378     aMAdd1.Clear();
379     aNbAdd = aMAdd.Extent();
380     for (i=1; i<=aNbAdd; ++i) {
381       const TopoDS_Shape& aF=aMAdd(i);
382       //
383       //aMAdd1.Clear();
384       aExp.Init(aF, TopAbs_EDGE);
385       for (; aExp.More(); aExp.Next()) {
386         const TopoDS_Shape& aE=aExp.Current();
387         if (theMEAvoid.Contains(aE)){
388           continue;
389         }
390         //
391         const BOPCol_ListOfShape& aLF=aMEF.FindFromKey(aE);
392         aIt.Initialize(aLF);
393         for (; aIt.More(); aIt.Next()) {
394           const TopoDS_Shape& aFx=aIt.Value();
395           if (aFx.IsSame(aF)) {
396             continue;
397           }
398           if (aMCB.Contains(aFx)) {
399             continue;
400           }
401           aMAdd1.Add(aFx);
402         }
403       }//for (; aExp.More(); aExp.Next()){
404       aMCB.Add(aF);
405     }// for (i=1; i<=aNbAdd; ++i) {
406     //
407     aNbAdd1=aMAdd1.Extent();
408     if (!aNbAdd1) {
409       break;
410     }
411     //
412     aMAdd.Clear();
413     for (i=1; i<=aNbAdd1; ++i) {
414       const TopoDS_Shape& aFAdd=aMAdd1(i);
415       aMAdd.Add(aFAdd);
416     }
417     //
418   }//while(1) {
419   
420   //
421   aNbF=aMCB.Extent();
422   for (i=1; i<=aNbF; ++i) {
423     const TopoDS_Shape& aF=aMCB(i);
424     theLCB.Append(aF);
425   }
426 }
427 //=======================================================================
428 // function:  ComputeStateByOnePoint
429 // purpose: 
430 //=======================================================================
431 TopAbs_State BOPTools_AlgoTools::ComputeStateByOnePoint(const TopoDS_Shape& theS,
432                                                         const TopoDS_Solid& theRef,
433                                                         const Standard_Real theTol,
434                                                         Handle(BOPInt_Context)& theContext)
435 {
436   TopAbs_State aState;
437   TopAbs_ShapeEnum aType;
438   //
439   aState=TopAbs_UNKNOWN;
440   aType=theS.ShapeType();
441   if (aType==TopAbs_VERTEX) {
442     const TopoDS_Vertex& aV=(*(TopoDS_Vertex*)(&theS));
443     aState=BOPTools_AlgoTools::ComputeState(aV, theRef, theTol, theContext);
444   }
445   else if (aType==TopAbs_EDGE) {
446     const TopoDS_Edge& aE=(*(TopoDS_Edge*)(&theS));
447     aState=BOPTools_AlgoTools::ComputeState(aE, theRef, theTol, theContext);
448   }
449   return aState;
450 }
451
452 //=======================================================================
453 // function:  ComputeState
454 // purpose: 
455 //=======================================================================
456 TopAbs_State BOPTools_AlgoTools::ComputeState(const TopoDS_Face& theF,
457                                               const TopoDS_Solid& theRef,
458                                               const Standard_Real theTol,
459                                               BOPCol_IndexedMapOfShape& theBounds,
460                                               Handle(BOPInt_Context)& theContext)
461 {
462   TopAbs_State aState;
463   TopExp_Explorer aExp; 
464   TopoDS_Edge aE1;
465   gp_Pnt2d aP2D;
466   gp_Pnt aP3D; 
467   //
468   aState=TopAbs_UNKNOWN;
469   //
470   aExp.Init(theF, TopAbs_EDGE);
471   for (; aExp.More(); aExp.Next()) {
472     const TopoDS_Edge& aSE=(*(TopoDS_Edge*)(&aExp.Current()));
473     if (BRep_Tool::Degenerated(aSE)) {
474       continue;
475     }
476     //
477     if (!theBounds.Contains(aSE)) {
478       const TopoDS_Edge& aE=(*(TopoDS_Edge*)(&aSE));
479       aState=BOPTools_AlgoTools::ComputeState(aE, theRef, theTol, theContext);
480       return aState;
481     }
482     if (aE1.IsNull()) {
483       aE1=(*(TopoDS_Edge*)(&aSE));
484     }
485   }
486   // !!<- process edges that are all on theRef
487   if (!aE1.IsNull()) {
488     BOPTools_AlgoTools3D::PointNearEdge(aE1, theF, aP2D, aP3D, theContext);
489     aState=BOPTools_AlgoTools::ComputeState(aP3D, theRef, theTol, theContext);
490   }
491   //
492   return aState;
493 }
494 //=======================================================================
495 // function:  ComputeState
496 // purpose: 
497 //=======================================================================
498   TopAbs_State BOPTools_AlgoTools::ComputeState(const TopoDS_Vertex& theV,
499                                             const TopoDS_Solid& theRef,
500                                             const Standard_Real theTol,
501                                             Handle(BOPInt_Context)& theContext)
502 {
503   TopAbs_State aState;
504   gp_Pnt aP3D; 
505   //
506   aP3D=BRep_Tool::Pnt(theV);
507   aState=BOPTools_AlgoTools::ComputeState(aP3D, theRef, theTol, theContext);
508   return aState;
509 }
510 //=======================================================================
511 // function:  ComputeState
512 // purpose: 
513 //=======================================================================
514   TopAbs_State BOPTools_AlgoTools::ComputeState(const TopoDS_Edge& theE,
515                                             const TopoDS_Solid& theRef,
516                                             const Standard_Real theTol,
517                                             Handle(BOPInt_Context)& theContext)
518 {
519   Standard_Real aT1, aT2, aT = 0.;
520   TopAbs_State aState;
521   Handle(Geom_Curve) aC3D;
522   gp_Pnt aP3D; 
523   //
524   aC3D = BRep_Tool::Curve(theE, aT1, aT2);
525   //
526   if(aC3D.IsNull()) {
527     //it means that we are in degenerated edge
528     const TopoDS_Vertex& aV = TopExp::FirstVertex(theE);
529     if(aV.IsNull()){
530       return TopAbs_UNKNOWN;
531     }
532     aP3D=BRep_Tool::Pnt(aV);
533   }
534   else {//usual case
535     Standard_Boolean bF2Inf, bL2Inf;
536     Standard_Real dT=10.;
537     //
538     bF2Inf = Precision::IsNegativeInfinite(aT1);
539     bL2Inf = Precision::IsPositiveInfinite(aT2);
540     //
541     if (bF2Inf && !bL2Inf) {
542       aT=aT2-dT;
543     }
544     else if (!bF2Inf && bL2Inf) {
545       aT=aT1+dT;
546     }
547     else if (bF2Inf && bL2Inf) {
548       aT=0.;
549     }
550     else {
551       aT=IntTools_Tools::IntermediatePoint(aT1, aT2);
552     }
553     aC3D->D0(aT, aP3D);
554   }
555   //
556   aState=BOPTools_AlgoTools::ComputeState(aP3D, theRef, theTol, theContext);
557   //
558   return aState;
559 }
560 //=======================================================================
561 // function:  ComputeState
562 // purpose: 
563 //=======================================================================
564   TopAbs_State BOPTools_AlgoTools::ComputeState(const gp_Pnt& theP,
565                                             const TopoDS_Solid& theRef,
566                                             const Standard_Real theTol,
567                                             Handle(BOPInt_Context)& theContext)
568 {
569   TopAbs_State aState;
570   //
571   BRepClass3d_SolidClassifier& aSC=theContext->SolidClassifier(theRef);
572   aSC.Perform(theP, theTol);
573   //
574   aState=aSC.State();
575   //
576   return aState;
577 }
578 //=======================================================================
579 //function : IsInternalFace
580 //purpose  : 
581 //=======================================================================
582 Standard_Integer BOPTools_AlgoTools::IsInternalFace
583   (const TopoDS_Face& theFace,
584    const TopoDS_Solid& theSolid,
585    BOPCol_IndexedDataMapOfShapeListOfShape& theMEF,
586    const Standard_Real theTol,
587    Handle(BOPInt_Context)& theContext)
588 {
589   Standard_Boolean bDegenerated;
590   Standard_Integer aNbF, iRet, iFound;
591   TopAbs_Orientation aOr;
592   TopoDS_Edge aE1;
593   TopExp_Explorer aExp;
594   BOPCol_ListIteratorOfListOfShape aItF;
595   //
596   // For all invoked functions: [::IsInternalFace(...)] 
597   // the returned value iRet means:
598   // iRet=0;  - state is not IN
599   // iRet=1;  - state is IN
600   // iRet=2;  - state can not be found by the method of angles
601   //
602   // For this function the returned value iRet means:
603   // iRet=0;  - state is not IN
604   // iRet=1;  - state is IN
605   //
606   iRet=0; 
607   // 1 Try to find an edge from theFace in theMEF
608   iFound=0;
609   aExp.Init(theFace, TopAbs_EDGE);
610   for(; aExp.More(); aExp.Next()) {
611     const TopoDS_Edge& aE=(*(TopoDS_Edge*)(&aExp.Current()));
612     if (!theMEF.Contains(aE)) {
613       continue;
614     }
615     //
616     ++iFound;
617     //
618     aOr=aE.Orientation();
619     if (aOr==TopAbs_INTERNAL) {
620       continue;
621     }
622     bDegenerated=BRep_Tool::Degenerated(aE);
623     if (bDegenerated){
624       continue;
625     }
626     // aE
627     BOPCol_ListOfShape& aLF=theMEF.ChangeFromKey(aE);
628     aNbF=aLF.Extent();
629     if (!aNbF) {
630       return iRet; // it can not be so
631     }
632     //
633     else if (aNbF==1) {
634       // aE is internal edge on aLF.First()
635       const TopoDS_Face& aF1=(*(TopoDS_Face*)(&aLF.First()));
636       BOPTools_AlgoTools::GetEdgeOnFace(aE, aF1, aE1);
637       if (aE1.Orientation()!=TopAbs_INTERNAL) {
638         iRet=2; 
639         break;
640       }
641       //
642       iRet=BOPTools_AlgoTools::IsInternalFace(theFace, aE, aF1, aF1, theContext);
643       break;
644     }
645     //
646     else if (aNbF==2) {
647       const TopoDS_Face& aF1=(*(TopoDS_Face*)(&aLF.First()));
648       const TopoDS_Face& aF2=(*(TopoDS_Face*)(&aLF.Last()));
649       //
650       if (aF2.IsSame(aF1) && BRep_Tool::IsClosed(aE, aF1)) {
651         // treat as it was for 1 face
652         iRet=BOPTools_AlgoTools::IsInternalFace(theFace, aE, aF1, aF2, theContext);
653         break;
654       }
655     }
656     //
657     if (aNbF%2) {
658       iRet=0;
659       return iRet; // it can not be so
660     }
661     else { // aNbF=2,4,6,8,...
662       iRet=BOPTools_AlgoTools::IsInternalFace(theFace, aE, aLF, theContext);
663       break;
664     }
665   }//for(; aExp.More(); aExp.Next()) {
666   //
667   if (!iFound) {
668     // the face has no shared edges with the solid
669     iRet=2;
670   }
671   //
672   if (iRet!=2) {
673     return iRet;
674   }
675   //
676   //========================================
677   // 2. Classify face using classifier
678   //
679   TopAbs_State aState;
680   BOPCol_IndexedMapOfShape aBounds;
681   //
682   BOPTools::MapShapes(theSolid, TopAbs_EDGE, aBounds);
683   //
684   aState=BOPTools_AlgoTools::ComputeState(theFace, theSolid, theTol, aBounds, theContext);
685   //
686   iRet=(aState==TopAbs_IN)? 1 : 0;
687   //
688   return iRet;
689 }
690 //=======================================================================
691 //function : IsInternalFace
692 //purpose  : 
693 //=======================================================================
694 Standard_Integer BOPTools_AlgoTools::IsInternalFace(const TopoDS_Face& theFace,
695                                                     const TopoDS_Edge& theEdge,
696                                                     BOPCol_ListOfShape& theLF,
697                                                     Handle(BOPInt_Context)& theContext)
698 {
699   Standard_Integer aNbF, iRet;
700   //
701   iRet=0;
702   //
703   aNbF=theLF.Extent();
704   if (aNbF==2) {
705     const TopoDS_Face& aF1=(*(TopoDS_Face*)(&theLF.First()));
706     const TopoDS_Face& aF2=(*(TopoDS_Face*)(&theLF.Last()));
707     iRet=BOPTools_AlgoTools::IsInternalFace(theFace, theEdge, aF1, aF2, theContext);
708     return iRet;
709   }
710   //
711   else {
712     BOPTools_ListOfCoupleOfShape aLCFF;
713     BOPTools_ListIteratorOfListOfCoupleOfShape aIt;
714     //
715     FindFacePairs(theEdge, theLF, aLCFF, theContext);
716     //
717     aIt.Initialize(aLCFF);
718     for (; aIt.More(); aIt.Next()) {
719       BOPTools_CoupleOfShape& aCSFF=aIt.ChangeValue();
720       //
721       const TopoDS_Face& aF1=(*(TopoDS_Face*)(&aCSFF.Shape1()));
722       const TopoDS_Face& aF2=(*(TopoDS_Face*)(&aCSFF.Shape2()));
723       iRet=BOPTools_AlgoTools::IsInternalFace(theFace, theEdge, aF1, aF2, theContext);
724       if (iRet) {
725         return iRet;
726       }
727     }
728   }
729   return iRet;
730 }
731 //=======================================================================
732 //function : IsInternalFace
733 //purpose  : 
734 //=======================================================================
735 Standard_Integer BOPTools_AlgoTools::IsInternalFace
736   (const TopoDS_Face& theFace,
737    const TopoDS_Edge& theEdge,
738    const TopoDS_Face& theFace1,
739    const TopoDS_Face& theFace2,
740    Handle(BOPInt_Context)& theContext)
741 {
742   Standard_Boolean bRet;
743   Standard_Integer iRet;
744   TopoDS_Edge aE1, aE2;
745   TopoDS_Face aFOff;
746   BOPTools_ListOfCoupleOfShape theLCSOff;
747   BOPTools_CoupleOfShape aCS1, aCS2;
748   //
749   BOPTools_AlgoTools::GetEdgeOnFace(theEdge, theFace1, aE1);
750   if (aE1.Orientation()==TopAbs_INTERNAL) {
751     aE2=aE1;
752     aE1.Orientation(TopAbs_FORWARD);
753     aE2.Orientation(TopAbs_REVERSED);
754   }
755   else if (theFace1==theFace2) {
756     aE2=aE1;
757     aE1.Orientation(TopAbs_FORWARD);
758     aE2.Orientation(TopAbs_REVERSED);
759   }
760   else {
761     BOPTools_AlgoTools::GetEdgeOnFace(theEdge, theFace2, aE2);
762   }
763   //
764   aCS1.SetShape1(theEdge);
765   aCS1.SetShape2(theFace);
766   theLCSOff.Append(aCS1);
767   //
768   aCS2.SetShape1(aE2);
769   aCS2.SetShape2(theFace2);
770   theLCSOff.Append(aCS2);
771   //
772   bRet=GetFaceOff(aE1, theFace1, theLCSOff, aFOff, theContext);
773   //
774   iRet=0; // theFace is not internal
775   if (theFace.IsEqual(aFOff)) {
776     // theFace is internal
777     iRet=1;  
778     if (!bRet) {
779       // theFace seems to be internal  
780       iRet=2;
781     }
782   }
783   return iRet;
784 }
785 //=======================================================================
786 //function : GetFaceOff
787 //purpose  : 
788 //=======================================================================
789 Standard_Boolean BOPTools_AlgoTools::GetFaceOff
790   (const TopoDS_Edge& theE1,
791    const TopoDS_Face& theF1,
792    BOPTools_ListOfCoupleOfShape& theLCSOff,
793    TopoDS_Face& theFOff,
794    Handle(BOPInt_Context)& theContext)
795 {
796   Standard_Boolean bRet;
797   Standard_Real aT, aT1, aT2, aAngle, aTwoPI, aAngleMin;
798   gp_Pnt aPn1, aPn2, aPx;
799   gp_Dir aDN1, aDN2, aDBF, aDBF2, aDTF;
800   gp_Vec aVTgt;
801   TopAbs_Orientation aOr;
802   Handle(Geom_Curve)aC3D;
803   BOPTools_ListIteratorOfListOfCoupleOfShape aIt;
804   //
805   aAngleMin=100.;
806   aTwoPI=M_PI+M_PI;
807   aC3D =BRep_Tool::Curve(theE1, aT1, aT2);
808   aT=BOPTools_AlgoTools2D::IntermediatePoint(aT1, aT2);
809   aC3D->D0(aT, aPx);
810   //
811   BOPTools_AlgoTools2D::EdgeTangent(theE1, aT, aVTgt);
812   gp_Dir aDTgt(aVTgt), aDTgt2;
813   aOr = theE1.Orientation();
814   //
815   GetFaceDir(theE1, theF1, aPx, aT, aDTgt, aDN1, aDBF, theContext);
816   //
817   aDTF=aDN1^aDBF;
818   //
819   bRet=Standard_True;
820   aIt.Initialize(theLCSOff);
821   for (; aIt.More(); aIt.Next()) {
822     const BOPTools_CoupleOfShape& aCS=aIt.Value();
823     const TopoDS_Edge& aE2=(*(TopoDS_Edge*)(&aCS.Shape1()));
824     const TopoDS_Face& aF2=(*(TopoDS_Face*)(&aCS.Shape2()));
825     //
826     aDTgt2 = (aE2.Orientation()==aOr) ? aDTgt : aDTgt.Reversed();
827     GetFaceDir(aE2, aF2, aPx, aT, aDTgt2, aDN2, aDBF2, theContext);
828     //Angle
829     aAngle=AngleWithRef(aDBF, aDBF2, aDTF);
830     //
831     if(aAngle<0.) {
832       aAngle=aTwoPI+aAngle;
833     }
834     //
835     if (aAngle<Precision::Angular()) {
836       if (aF2==theF1) {
837         aAngle=M_PI;
838       }
839       else if (aF2.IsSame(theF1)) {
840         aAngle=aTwoPI;
841       }
842     }
843     //
844     if (aAngle<aAngleMin){
845       aAngleMin=aAngle;
846       theFOff=aF2;
847     }
848     else if (aAngle==aAngleMin) {
849       // the minimal angle can not be found
850       bRet=Standard_False;  
851     }
852   }
853   return bRet;
854 }
855 //=======================================================================
856 //function : GetEdgeOff
857 //purpose  : 
858 //=======================================================================
859   Standard_Boolean BOPTools_AlgoTools::GetEdgeOff(const TopoDS_Edge& theE1,
860                                               const TopoDS_Face& theF2,
861                                               TopoDS_Edge& theE2)
862 {
863   Standard_Boolean bFound;
864   TopAbs_Orientation aOr1, aOr1C, aOr2;
865   TopExp_Explorer anExp;
866   //
867   bFound=Standard_False;
868   aOr1=theE1.Orientation();
869   aOr1C=TopAbs::Reverse(aOr1);
870   //
871   anExp.Init(theF2, TopAbs_EDGE);
872   for (; anExp.More(); anExp.Next()) {
873     const TopoDS_Edge& aEF2=(*(TopoDS_Edge*)(&anExp.Current()));
874     if (aEF2.IsSame(theE1)) {
875       aOr2=aEF2.Orientation();
876       if (aOr2==aOr1C) {
877         theE2=aEF2;
878         bFound=!bFound;
879         return bFound;
880       }
881     }
882   }
883   return bFound;
884 }
885
886 //=======================================================================
887 //function : AreFacesSameDomain
888 //purpose  : 
889 //=======================================================================
890   Standard_Boolean BOPTools_AlgoTools::AreFacesSameDomain(const TopoDS_Face& theF1,
891                                                           const TopoDS_Face& theF2,
892                                                           Handle(BOPInt_Context)& theContext)
893 {
894   Standard_Boolean bFlag;
895   Standard_Integer iErr;
896   Standard_Real aTolF1, aTolF2, aTol;
897   gp_Pnt2d aP2D;
898   gp_Pnt aP;
899   TopoDS_Face aF1, aF2;
900   TopoDS_Edge aE1;
901   TopExp_Explorer aExp;
902   //
903   bFlag=Standard_False;
904   //
905   aF1=theF1;
906   aF1.Orientation(TopAbs_FORWARD);
907   aF2=theF2;
908   aF2.Orientation(TopAbs_FORWARD);
909   //
910   aTolF1=BRep_Tool::Tolerance(aF1);
911   // 1
912   aExp.Init(aF1, TopAbs_EDGE);
913   for (; aExp.More(); aExp.Next()) {
914     aE1=(*(TopoDS_Edge*)(&aExp.Current()));
915     if (!BRep_Tool::Degenerated(aE1)) {
916       Standard_Real aTolE = BRep_Tool::Tolerance(aE1);
917       aTolF1 = (aTolE > aTolF1) ? aTolE : aTolF1;
918     }
919   }
920   // 2
921   aTolF2=BRep_Tool::Tolerance(aF2);
922   aTol=aTolF1+aTolF2;
923   //
924   iErr = BOPTools_AlgoTools3D::PointInFace(aF1, aP, aP2D, theContext);
925   if (!iErr) {
926     bFlag=theContext->IsValidPointForFace(aP, aF2, aTol);
927   }
928   //
929   return bFlag;
930 }
931
932 //=======================================================================
933 //function : CheckSameGeom
934 //purpose  : 
935 //=======================================================================
936   Standard_Boolean BOPTools_AlgoTools::CheckSameGeom(const TopoDS_Face& theF1,
937                                                  const TopoDS_Face& theF2,
938                                                  Handle(BOPInt_Context)& theContext)
939 {
940   Standard_Boolean bRet;
941   Standard_Real aTolF1, aTolF2, aTol;
942   gp_Pnt2d aP2D;
943   gp_Pnt aP;
944   TopExp_Explorer aExp;
945   //
946   bRet=Standard_False;
947   aExp.Init(theF1, TopAbs_EDGE);
948   for (; aExp.More(); aExp.Next()) {
949     const TopoDS_Edge& aE=(*(TopoDS_Edge*)(&aExp.Current()));
950     if (!BRep_Tool::Degenerated(aE)) {
951       aTolF1=BRep_Tool::Tolerance(theF1);
952       aTolF2=BRep_Tool::Tolerance(theF2);
953       aTol=aTolF1+aTolF2;
954       BOPTools_AlgoTools3D::PointNearEdge(aE, theF1, aP2D, aP, theContext);
955       bRet=theContext->IsValidPointForFace(aP, theF2, aTol);
956       break;
957     }
958   }
959   return bRet;
960 }
961 //=======================================================================
962 // function: Sense
963 // purpose: 
964 //=======================================================================
965   Standard_Integer BOPTools_AlgoTools::Sense (const TopoDS_Face& theF1,
966                                           const TopoDS_Face& theF2)
967 {
968   Standard_Integer iSense=0;
969   gp_Dir aDNF1, aDNF2;
970   TopoDS_Edge aE1, aE2;
971   TopExp_Explorer aExp;
972   //
973   aExp.Init(theF1, TopAbs_EDGE);
974   for (; aExp.More(); aExp.Next()) {
975     aE1=(*(TopoDS_Edge*)(&aExp.Current()));
976     if (!BRep_Tool::Degenerated(aE1)) {
977       if (!BRep_Tool::IsClosed(aE1, theF1)) {
978         break;
979       }
980     }
981   }
982   //
983   aExp.Init(theF2, TopAbs_EDGE);
984   for (; aExp.More(); aExp.Next()) {
985     aE2=(*(TopoDS_Edge*)(&aExp.Current()));
986     if (!BRep_Tool::Degenerated(aE2)) {
987       if (!BRep_Tool::IsClosed(aE2, theF2)) {
988         if (aE2.IsSame(aE1)) {
989           iSense=1;
990           break;
991         }
992       }
993     }
994   }
995   //
996   if (!iSense) {
997     return iSense;
998   }
999   //
1000   BOPTools_AlgoTools3D::GetNormalToFaceOnEdge(aE1, theF1, aDNF1);
1001   BOPTools_AlgoTools3D::GetNormalToFaceOnEdge(aE2, theF2, aDNF2);
1002   //
1003   iSense=BOPTools_AlgoTools3D::SenseFlag(aDNF1, aDNF2);
1004   //
1005   return iSense;
1006 }
1007 //=======================================================================
1008 // function: IsSplitToReverse
1009 // purpose: 
1010 //=======================================================================
1011   Standard_Boolean BOPTools_AlgoTools::IsSplitToReverse(const TopoDS_Shape& theSp,
1012                                                     const TopoDS_Shape& theSr,
1013                                                     Handle(BOPInt_Context)& theContext)
1014 {
1015   Standard_Boolean bRet;
1016   TopAbs_ShapeEnum aType;
1017   //
1018   bRet=Standard_False;
1019   //
1020   aType=theSp.ShapeType();
1021   switch (aType) {
1022     case TopAbs_EDGE: {
1023       const TopoDS_Edge& aESp=(*(TopoDS_Edge*)(&theSp));
1024       const TopoDS_Edge& aESr=(*(TopoDS_Edge*)(&theSr));
1025       bRet=BOPTools_AlgoTools::IsSplitToReverse(aESp, aESr, theContext);
1026     }
1027       break;
1028       //
1029     case TopAbs_FACE: {
1030       const TopoDS_Face& aFSp=(*(TopoDS_Face*)(&theSp));
1031       const TopoDS_Face& aFSr=(*(TopoDS_Face*)(&theSr));
1032       bRet=BOPTools_AlgoTools::IsSplitToReverse(aFSp, aFSr, theContext);
1033     }
1034       break;
1035       //
1036     default:
1037       break;
1038   }
1039   return bRet;
1040 }
1041 //=======================================================================
1042 //function :IsSplitToReverse
1043 //purpose  : 
1044 //=======================================================================
1045   Standard_Boolean BOPTools_AlgoTools::IsSplitToReverse(const TopoDS_Face& theFSp,
1046                                                     const TopoDS_Face& theFSr,
1047                                                     Handle(BOPInt_Context)& theContext)
1048 {
1049   Standard_Boolean bRet, bFound, bInFace;
1050   Standard_Real aT1, aT2, aT, aU, aV, aScPr;
1051   gp_Pnt aPFSp, aPFSr;
1052   gp_Dir aDNFSp;
1053   gp_Vec aD1U, aD1V;
1054   Handle(Geom_Surface) aSr, aSp;
1055   TopAbs_Orientation aOrSr, aOrSp;
1056   TopExp_Explorer anExp;
1057   TopoDS_Edge aESp;
1058   //
1059   bRet=Standard_False;
1060   //
1061   aSr=BRep_Tool::Surface(theFSr);
1062   aSp=BRep_Tool::Surface(theFSp);
1063   if (aSr==aSp) {
1064     aOrSr=theFSr.Orientation();
1065     aOrSp=theFSp.Orientation();
1066     bRet=(aOrSr!=aOrSp);
1067     return bRet;
1068   }
1069   //
1070   bFound=Standard_False;
1071   anExp.Init(theFSp, TopAbs_EDGE);
1072   for (; anExp.More(); anExp.Next()) {
1073     aESp=(*(TopoDS_Edge*)(&anExp.Current()));
1074     if (!BRep_Tool::Degenerated(aESp)) {
1075       if (!BRep_Tool::IsClosed(aESp, theFSp)) {
1076         bFound=!bFound;
1077         break;
1078       }
1079     }
1080   }
1081   if (!bFound) {
1082     Standard_Boolean bFlag;
1083     Standard_Integer iErr;
1084     gp_Pnt2d aP2DFSp;
1085     //
1086     iErr=BOPTools_AlgoTools3D::PointInFace(theFSp, aPFSp, aP2DFSp, theContext);
1087     if (iErr) {
1088       return bRet;
1089     }
1090     //
1091     aP2DFSp.Coord(aU, aV);
1092     bFlag=BOPTools_AlgoTools3D::GetNormalToSurface(aSp, aU, aV, aDNFSp);
1093     if (!bFlag) {
1094       return bRet;
1095     }
1096   }
1097   else {
1098     BRep_Tool::Range(aESp, aT1, aT2);
1099     aT=BOPTools_AlgoTools2D::IntermediatePoint(aT1, aT2);
1100     BOPTools_AlgoTools3D::GetApproxNormalToFaceOnEdge(aESp, theFSp, aT, aPFSp, aDNFSp, theContext);
1101   }
1102   //
1103   // Parts of theContext->ComputeVS(..) 
1104   GeomAPI_ProjectPointOnSurf& aProjector=theContext->ProjPS(theFSr);
1105   aProjector.Perform(aPFSp);
1106   if (!aProjector.IsDone()) {
1107     return bRet;
1108   }
1109   //
1110   aProjector.LowerDistanceParameters(aU, aV);
1111   gp_Pnt2d aP2D(aU, aV);
1112   bInFace=theContext->IsPointInFace (theFSr, aP2D);
1113   if (!bInFace) {
1114     return bRet;
1115   }
1116   //
1117   aSr->D1(aU, aV, aPFSr, aD1U, aD1V);
1118   gp_Dir aDD1U(aD1U); 
1119   gp_Dir aDD1V(aD1V);
1120   gp_Dir aDNFSr=aDD1U^aDD1V; 
1121   if (theFSr.Orientation()==TopAbs_REVERSED){
1122     aDNFSr.Reverse();
1123   }
1124   //
1125   aScPr=aDNFSp*aDNFSr;
1126   bRet=(aScPr<0.);
1127   //
1128   return bRet;
1129 }
1130 //=======================================================================
1131 //function :IsSplitToReverse
1132 //purpose  : 
1133 //=======================================================================
1134   Standard_Boolean BOPTools_AlgoTools::IsSplitToReverse(const TopoDS_Edge& aEF1,
1135                                                     const TopoDS_Edge& aEF2,
1136                                                     Handle(BOPInt_Context)& theContext)
1137 {
1138   Standard_Boolean bRet, bIsDegenerated;
1139   //
1140   bRet=Standard_False;
1141   bIsDegenerated=(BRep_Tool::Degenerated(aEF1) || 
1142                   BRep_Tool::Degenerated(aEF2));
1143   if (bIsDegenerated) {
1144     return bRet;
1145   }
1146   //
1147   Standard_Real a, b;
1148   TopAbs_Orientation aOrE, aOrSp;
1149   Handle(Geom_Curve)aC1, aC2;
1150   //
1151   aC2=BRep_Tool::Curve(aEF2, a, b);
1152   aC1=BRep_Tool::Curve(aEF1, a, b);
1153   //
1154   if (aC1==aC2) {
1155     aOrE=aEF2.Orientation();
1156     aOrSp=aEF1.Orientation();
1157     bRet=(aOrE!=aOrSp);
1158     return bRet;
1159   }
1160   //
1161   Standard_Real aT1, aT2, aScPr;
1162   gp_Vec aV1, aV2;
1163   gp_Pnt aP;
1164   //
1165   aT1=BOPTools_AlgoTools2D::IntermediatePoint(a, b);
1166   aC1->D0(aT1, aP);
1167   BOPTools_AlgoTools2D::EdgeTangent(aEF1, aT1, aV1);
1168   gp_Dir aDT1(aV1);
1169   //
1170   theContext->ProjectPointOnEdge(aP, aEF2, aT2);
1171   //
1172   BOPTools_AlgoTools2D::EdgeTangent(aEF2, aT2, aV2);
1173   gp_Dir aDT2(aV2);
1174   //
1175   aScPr=aDT1*aDT2;
1176   bRet=(aScPr<0.);
1177   //
1178   return bRet;
1179 }
1180
1181 //=======================================================================
1182 //function : IsHole
1183 //purpose  : 
1184 //=======================================================================
1185   Standard_Boolean BOPTools_AlgoTools::IsHole(const TopoDS_Shape& aW,
1186                                           const TopoDS_Shape& aFace)
1187 {
1188   Standard_Boolean bIsHole;
1189   Standard_Integer i, aNbS;
1190   Standard_Real aT1, aT2, aS;
1191   Standard_Real aU1, aU2, aU, dU;
1192   Standard_Real aX1, aY1, aX0, aY0;
1193   TopAbs_Orientation aOr;
1194   
1195   gp_Pnt2d aP2D0, aP2D1;
1196   Handle(Geom2d_Curve) aC2D;
1197   TopoDS_Face aF, aFF;
1198   TopoDS_Iterator aItW;
1199   //
1200   bIsHole=Standard_False;
1201   //
1202   aF=(*(TopoDS_Face *)(&aFace));
1203   aFF=aF;
1204   aFF.Orientation(TopAbs_FORWARD);
1205   //
1206   aS=0.;
1207   aItW.Initialize(aW);
1208   for (; aItW.More(); aItW.Next()) { 
1209     const TopoDS_Edge& aE=(*(TopoDS_Edge *)(&aItW.Value()));
1210     aOr=aE.Orientation();
1211     if (!(aOr==TopAbs_FORWARD || 
1212           aOr==TopAbs_REVERSED)) {
1213       continue;
1214     }
1215     //
1216     aC2D=BRep_Tool::CurveOnSurface(aE, aFF, aT1, aT2);
1217     if (aC2D.IsNull()) {
1218       break; //xx
1219     }
1220     //
1221     BRepAdaptor_Curve2d aBAC2D(aE, aFF);
1222     aNbS=Geom2dInt_Geom2dCurveTool::NbSamples(aBAC2D);
1223     if (aNbS>2) {
1224       aNbS*=4;
1225     }
1226     //
1227     dU=(aT2-aT1)/(Standard_Real)(aNbS-1);
1228     aU =aT1;
1229     aU1=aT1;
1230     aU2=aT2;
1231     if (aOr==TopAbs_REVERSED) {
1232       aU =aT2;
1233       aU1=aT2;
1234       aU2=aT1;
1235       dU=-dU;
1236     }
1237     //
1238     aC2D->D0(aU, aP2D0);
1239     for(i=2; i<=aNbS; i++) {
1240       aU=aU1+(i-1)*dU;
1241       aC2D->D0(aU, aP2D1);
1242       aP2D0.Coord(aX0, aY0);
1243       aP2D1.Coord(aX1, aY1);
1244       //
1245       aS=aS+(aY0+aY1)*(aX1-aX0); 
1246       //
1247       aP2D0=aP2D1;
1248     }
1249   }//for (; aItW.More(); aItW.Next()) { 
1250   bIsHole=(aS>0.);
1251   return bIsHole;
1252 }
1253
1254 //=======================================================================
1255 // function: MakeContainer
1256 // purpose: 
1257 //=======================================================================
1258   void BOPTools_AlgoTools::MakeContainer(const TopAbs_ShapeEnum theType,
1259                                      TopoDS_Shape& theC)
1260 {
1261   BRep_Builder aBB;
1262   //
1263   switch(theType) {
1264     case TopAbs_COMPOUND:{
1265       TopoDS_Compound aC;
1266       aBB.MakeCompound(aC);
1267       theC=aC;
1268     }
1269       break;
1270       //
1271     case TopAbs_COMPSOLID:{
1272       TopoDS_CompSolid aCS;
1273       aBB.MakeCompSolid(aCS);
1274       theC=aCS;
1275     }
1276       break;
1277       //
1278     case TopAbs_SOLID:{
1279       TopoDS_Solid aSolid;
1280       aBB.MakeSolid(aSolid);
1281       theC=aSolid;
1282     }  
1283       break;
1284       //
1285       //
1286     case TopAbs_SHELL:{
1287       TopoDS_Shell aShell;
1288       aBB.MakeShell(aShell);
1289       theC=aShell;
1290     }  
1291       break;
1292       //
1293     case TopAbs_WIRE: {
1294       TopoDS_Wire aWire;
1295       aBB.MakeWire(aWire);
1296       theC=aWire;
1297     }
1298       break;
1299       //
1300     default:
1301       break;
1302   }
1303 }
1304 //=======================================================================
1305 // function: MakePCurve
1306 // purpose: 
1307 //=======================================================================
1308   void  BOPTools_AlgoTools::MakePCurve(const TopoDS_Edge& aE,
1309                                        const TopoDS_Face& aF1,
1310                                        const TopoDS_Face& aF2,
1311                                        const IntTools_Curve& aIC,
1312                                        const Standard_Boolean bPC1,
1313                                        const Standard_Boolean bPC2)
1314
1315 {
1316   Standard_Integer i;
1317   Standard_Real aTolE, aT1, aT2, aOutFirst, aOutLast, aOutTol;
1318   Handle(Geom2d_Curve) aC2D, aC2DA, aC2Dx1;
1319   TopoDS_Face aFFWD; 
1320   BRep_Builder aBB;
1321   Standard_Boolean bPC;
1322   //
1323   aTolE=BRep_Tool::Tolerance(aE);
1324   //
1325   const Handle(Geom_Curve)& aC3DE=BRep_Tool::Curve(aE, aT1, aT2);
1326   Handle(Geom_TrimmedCurve)aC3DETrim=new Geom_TrimmedCurve(aC3DE, aT1, aT2);
1327   //
1328   for (i=0; i<2; ++i) {
1329     bPC = !i ? bPC1 : bPC2;
1330     if (!bPC) {
1331       continue;
1332     }
1333     //
1334     if (!i) {
1335       aFFWD=aF1;
1336       aC2Dx1=aIC.FirstCurve2d();
1337     }
1338     else {
1339       aFFWD=aF2;
1340       aC2Dx1=aIC.SecondCurve2d();
1341     }
1342     //
1343     aFFWD.Orientation(TopAbs_FORWARD);
1344     //
1345     aC2D=aC2Dx1;
1346     if (aC2D.IsNull()) { 
1347       BOPTools_AlgoTools2D::BuildPCurveForEdgeOnFace(aE, aFFWD);
1348       BOPTools_AlgoTools2D::CurveOnSurface(aE, aFFWD, aC2D, 
1349                                        aOutFirst, aOutLast, 
1350                                        aOutTol);
1351       }
1352     //
1353     if (aC3DE->IsPeriodic()) {
1354       BOPTools_AlgoTools2D::AdjustPCurveOnFace(aFFWD, aT1, aT2,  aC2D, aC2DA); 
1355     }
1356     else {
1357       BOPTools_AlgoTools2D::AdjustPCurveOnFace(aFFWD, aC3DETrim, aC2D, aC2DA); 
1358     }
1359     //
1360     aBB.UpdateEdge(aE, aC2DA, aFFWD, aTolE);
1361     //BRepLib::SameParameter(aE);
1362   }
1363   BRepLib::SameParameter(aE);
1364 }
1365 //=======================================================================
1366 // function: MakeEdge
1367 // purpose: 
1368 //=======================================================================
1369   void  BOPTools_AlgoTools::MakeEdge(const IntTools_Curve& theIC,
1370                                  const TopoDS_Vertex& theV1,
1371                                  const Standard_Real theT1,
1372                                  const TopoDS_Vertex& theV2,
1373                                  const Standard_Real theT2,
1374                                  const Standard_Real theTolR3D,
1375                                  TopoDS_Edge& theE)
1376 {
1377   Standard_Real aTolV;
1378   BRep_Builder aBB;
1379   //
1380   BOPTools_AlgoTools::MakeSectEdge (theIC, theV1, theT1, theV2, theT2, theE);
1381   //
1382   aBB.UpdateEdge(theE, theTolR3D);
1383   //
1384   aTolV=BRep_Tool::Tolerance(theV1);
1385   if (aTolV<theTolR3D) {
1386     aBB.UpdateVertex(theV1, theTolR3D);
1387   }
1388   //
1389   aTolV=BRep_Tool::Tolerance(theV2);
1390   if (aTolV<theTolR3D) {
1391     aBB.UpdateVertex(theV2, theTolR3D);
1392   }
1393 }
1394 //=======================================================================
1395 // function: ComputeVV
1396 // purpose: 
1397 //=======================================================================
1398   Standard_Integer BOPTools_AlgoTools::ComputeVV(const TopoDS_Vertex& aV1, 
1399                                              const gp_Pnt& aP2,
1400                                              const Standard_Real aTolP2)
1401 {
1402   Standard_Real aTolV1, aTolSum, aTolSum2, aD2;
1403   gp_Pnt aP1;
1404   //
1405   aTolV1=BRep_Tool::Tolerance(aV1);
1406   
1407   aTolSum=aTolV1+aTolP2;
1408   aTolSum2=aTolSum*aTolSum;
1409   //
1410   aP1=BRep_Tool::Pnt(aV1);
1411   //
1412   aD2=aP1.SquareDistance(aP2);
1413   if (aD2>aTolSum2) {
1414     return 1;
1415   }
1416   return 0;
1417
1418 //=======================================================================
1419 // function: ComputeVV
1420 // purpose: 
1421 //=======================================================================
1422   Standard_Integer BOPTools_AlgoTools::ComputeVV(const TopoDS_Vertex& aV1, 
1423                                              const TopoDS_Vertex& aV2)
1424 {
1425   Standard_Real aTolV1, aTolV2, aTolSum, aTolSum2, aD2;
1426   gp_Pnt aP1, aP2;
1427   //
1428   aTolV1=BRep_Tool::Tolerance(aV1);
1429   aTolV2=BRep_Tool::Tolerance(aV2);
1430   aTolSum=aTolV1+aTolV2;
1431   aTolSum2=aTolSum*aTolSum;
1432   //
1433   aP1=BRep_Tool::Pnt(aV1);
1434   aP2=BRep_Tool::Pnt(aV2);
1435   //
1436   aD2=aP1.SquareDistance(aP2);
1437   if (aD2>aTolSum2) {
1438     return 1;
1439   }
1440   return 0;
1441 }
1442 //=======================================================================
1443 // function: MakeVertex
1444 // purpose : 
1445 //=======================================================================
1446   void BOPTools_AlgoTools::MakeVertex(BOPCol_ListOfShape& aLV,
1447                                   TopoDS_Vertex& aVnew)
1448 {
1449   Standard_Integer aNb;
1450   Standard_Real aTi, aDi, aDmax;
1451   gp_Pnt aPi, aP;
1452   gp_XYZ aXYZ(0.,0.,0.), aXYZi;
1453   BOPCol_ListIteratorOfListOfShape aIt;
1454   //
1455   aNb=aLV.Extent();
1456   if (aNb) {
1457     aIt.Initialize(aLV);
1458     for (; aIt.More(); aIt.Next()) {
1459       TopoDS_Vertex& aVi=*((TopoDS_Vertex*)(&aIt.Value()));
1460       aPi=BRep_Tool::Pnt(aVi);
1461       aXYZi=aPi.XYZ();
1462       aXYZ=aXYZ+aXYZi;
1463     }
1464     //
1465     aXYZ.Divide((Standard_Real)aNb);
1466     aP.SetXYZ(aXYZ);
1467     //
1468     aDmax=-1.;
1469     aIt.Initialize(aLV);
1470     for (; aIt.More(); aIt.Next()) {
1471       TopoDS_Vertex& aVi=*((TopoDS_Vertex*)(&aIt.Value()));
1472       aPi=BRep_Tool::Pnt(aVi);
1473       aTi=BRep_Tool::Tolerance(aVi);
1474       aDi=aP.SquareDistance(aPi);
1475       aDi=fsqrt(aDi);
1476       aDi=aDi+aTi;
1477       if (aDi > aDmax) {
1478         aDmax=aDi;
1479       }
1480     }
1481     //
1482     BRep_Builder aBB;
1483     aBB.MakeVertex (aVnew, aP, aDmax);
1484   }
1485 }
1486 //=======================================================================
1487 //function : GetEdgeOnFace
1488 //purpose  : 
1489 //=======================================================================
1490   Standard_Boolean BOPTools_AlgoTools::GetEdgeOnFace(const TopoDS_Edge& theE1,
1491                                                  const TopoDS_Face& theF2,
1492                                                  TopoDS_Edge& theE2)
1493 {
1494   Standard_Boolean bFound;
1495   TopoDS_Iterator aItF, aItW;
1496   //
1497   bFound=Standard_False;
1498   //
1499   aItF.Initialize(theF2);
1500   for (; aItF.More(); aItF.Next()) {
1501     const TopoDS_Shape& aW=aItF.Value();
1502     aItW.Initialize(aW);
1503     for (; aItW.More(); aItW.Next()) {
1504       const TopoDS_Shape& aE=aItW.Value();
1505       if (aE.IsSame(theE1)) {
1506         theE2=(*(TopoDS_Edge*)(&aE));
1507         bFound=!bFound;
1508         return bFound;
1509       }
1510     }
1511   }
1512   return bFound;
1513 }
1514 //=======================================================================
1515 //function : FindFacePairs
1516 //purpose  : 
1517 //=======================================================================
1518 Standard_Boolean FindFacePairs (const TopoDS_Edge& theE,
1519                                 const BOPCol_ListOfShape& thLF,
1520                                 BOPTools_ListOfCoupleOfShape& theLCFF,
1521                                 Handle(BOPInt_Context)& theContext)
1522 {
1523   Standard_Boolean bFound;
1524   Standard_Integer i, aNbCEF;
1525   TopAbs_Orientation aOr, aOrC = TopAbs_FORWARD;
1526   BOPCol_MapOfShape aMFP;
1527   TopoDS_Face aF1, aF2;
1528   TopoDS_Edge aEL, aE1;
1529   BOPCol_ListIteratorOfListOfShape aItLF;
1530   BOPTools_CoupleOfShape aCEF, aCFF;
1531   BOPTools_ListOfCoupleOfShape aLCEF, aLCEFx;
1532   BOPTools_ListIteratorOfListOfCoupleOfShape aIt;
1533   //
1534   bFound=Standard_True;
1535   //
1536   // Preface aLCEF
1537   aItLF.Initialize(thLF);
1538   for (; aItLF.More(); aItLF.Next()) { 
1539     const TopoDS_Face& aFL=(*(TopoDS_Face*)(&aItLF.Value()));
1540     //
1541     bFound=BOPTools_AlgoTools::GetEdgeOnFace(theE, aFL, aEL);
1542     if (!bFound) {
1543       return bFound; // it can not be so
1544     }
1545     //
1546     aCEF.SetShape1(aEL);
1547     aCEF.SetShape2(aFL);
1548     aLCEF.Append(aCEF);
1549   }
1550   //
1551   aNbCEF=aLCEF.Extent();
1552   while(aNbCEF) {
1553     //
1554     // aLCEFx
1555     aLCEFx.Clear();
1556     aIt.Initialize(aLCEF);
1557     for (i=0; aIt.More(); aIt.Next(), ++i) {
1558       const BOPTools_CoupleOfShape& aCSx=aIt.Value();
1559       const TopoDS_Shape& aEx=aCSx.Shape1();
1560       const TopoDS_Shape& aFx=aCSx.Shape2();
1561       //
1562       aOr=aEx.Orientation();
1563       //
1564       if (!i) {
1565         aOrC=TopAbs::Reverse(aOr);
1566         aE1=(*(TopoDS_Edge*)(&aEx));
1567         aF1=(*(TopoDS_Face*)(&aFx));
1568         aMFP.Add(aFx);
1569         continue;
1570       }
1571       //
1572       if (aOr==aOrC) {
1573         aLCEFx.Append(aCSx);
1574         aMFP.Add(aFx);
1575       }
1576     }
1577     //
1578     // F2
1579     BOPTools_AlgoTools::GetFaceOff(aE1, aF1, aLCEFx, aF2, theContext);
1580     //
1581     aCFF.SetShape1(aF1);
1582     aCFF.SetShape2(aF2);
1583     theLCFF.Append(aCFF);
1584     //
1585     aMFP.Add(aF1);
1586     aMFP.Add(aF2);
1587     //
1588     // refine aLCEF
1589     aLCEFx.Clear();
1590     aLCEFx=aLCEF;
1591     aLCEF.Clear();
1592     aIt.Initialize(aLCEFx);
1593     for (; aIt.More(); aIt.Next()) {
1594       const BOPTools_CoupleOfShape& aCSx=aIt.Value();
1595       const TopoDS_Shape& aFx=aCSx.Shape2();
1596       if (!aMFP.Contains(aFx)) {
1597         aLCEF.Append(aCSx);
1598       }
1599     }
1600     //
1601     aNbCEF=aLCEF.Extent();
1602   }//while(aNbCEF) {
1603   //
1604   return bFound;
1605 }
1606 //=======================================================================
1607 //function : AngleWithRef
1608 //purpose  : 
1609 //=======================================================================
1610 Standard_Real AngleWithRef(const gp_Dir& theD1,
1611                            const gp_Dir& theD2,
1612                            const gp_Dir& theDRef)
1613 {
1614   Standard_Real aCosinus, aSinus, aBeta, aHalfPI, aScPr;
1615   gp_XYZ aXYZ;
1616   //
1617   aHalfPI=0.5*M_PI;
1618   //
1619   const gp_XYZ& aXYZ1=theD1.XYZ();
1620   const gp_XYZ& aXYZ2=theD2.XYZ();
1621   aXYZ=aXYZ1.Crossed(aXYZ2);
1622   aSinus=aXYZ.Modulus();
1623   aCosinus=theD1*theD2;
1624   //
1625   aBeta=0.;
1626   if (aSinus>=0.) {
1627     aBeta=aHalfPI*(1.-aCosinus);
1628   }
1629   else {
1630     aBeta=2.*M_PI-aHalfPI*(3.+aCosinus);
1631   }
1632   //
1633   aScPr=aXYZ.Dot(theDRef.XYZ());
1634   if (aScPr<0.) {
1635     aBeta=-aBeta;
1636   }
1637   return aBeta;
1638 }
1639 //=======================================================================
1640 //function : fsqrt
1641 //purpose  : 
1642 //=======================================================================
1643 Standard_Real fsqrt(Standard_Real val)  
1644 {
1645   union {
1646     int tmp;
1647     float val;
1648   } u;
1649   //
1650   u.val = (float)val;
1651   u.tmp -= 1<<23; /* Remove last bit so 1.0 gives 1.0 */
1652   /* tmp is now an approximation to logbase2(val) */
1653   u.tmp >>= 1; /* divide by 2 */
1654   u.tmp += 1<<29; /* add 64 to exponent: (e+127)/2 =(e/2)+63, */
1655   /* that represents (e/2)-64 but we want e/2 */
1656   return (double)u.val;
1657 }
1658
1659 //=======================================================================
1660 // function: IsBlockInOnFace
1661 // purpose: 
1662 //=======================================================================
1663   Standard_Boolean BOPTools_AlgoTools::IsBlockInOnFace (const IntTools_Range& aShrR,
1664                                                         const TopoDS_Face& aF,
1665                                                         const TopoDS_Edge& aE1,
1666                                                         Handle(BOPInt_Context)& aContext)
1667 {
1668   Standard_Boolean bFlag;
1669   Standard_Real f1, l1, ULD, VLD;
1670   gp_Pnt2d aP2D;
1671   gp_Pnt aP11, aP12;
1672   //
1673   aShrR.Range(f1, l1);
1674   Standard_Real dt=0.0075,  k;//dt=0.001,  k;
1675   k=dt*(l1-f1);
1676   f1=f1+k;
1677   l1=l1-k;
1678   //
1679   // Treatment P11
1680   BOPTools_AlgoTools::PointOnEdge(aE1, f1, aP11);
1681   //
1682   GeomAPI_ProjectPointOnSurf& aProjector=aContext->ProjPS(aF);
1683   aProjector.Perform(aP11);
1684   //
1685   bFlag=aProjector.IsDone();
1686   if (!bFlag) {
1687     return bFlag;
1688   }
1689   
1690   aProjector.LowerDistanceParameters(ULD, VLD);
1691   aP2D.SetCoord(ULD, VLD);
1692   //
1693   bFlag=aContext->IsPointInOnFace (aF, aP2D);
1694   //
1695   if (!bFlag) {
1696     return bFlag;
1697   }
1698   //
1699   // Treatment P12
1700   BOPTools_AlgoTools::PointOnEdge(aE1, l1, aP12);
1701   //
1702   aProjector.Perform(aP12);
1703   //
1704   bFlag=aProjector.IsDone();
1705   if (!bFlag) {
1706     return bFlag;
1707   }
1708   
1709   aProjector.LowerDistanceParameters(ULD, VLD);
1710   aP2D.SetCoord(ULD, VLD);
1711   //
1712   bFlag=aContext->IsPointInOnFace (aF, aP2D);
1713   //
1714   if (!bFlag) {
1715     return bFlag;
1716   }
1717   //
1718   // Treatment intemediate
1719   Standard_Real m1, aTolF, aTolE, aTol, aDist;
1720   m1=IntTools_Tools::IntermediatePoint(f1, l1);
1721   BOPTools_AlgoTools::PointOnEdge(aE1, m1, aP12);
1722   //
1723   aProjector.Perform(aP12);
1724   //
1725   bFlag=aProjector.IsDone();
1726   if (!bFlag) {
1727     return bFlag;
1728   }
1729   //
1730   aTolE=BRep_Tool::Tolerance(aE1);
1731   aTolF=BRep_Tool::Tolerance(aF);
1732   aTol=aTolE+aTolF;
1733   aDist=aProjector.LowerDistance();
1734   if (aDist > aTol){
1735    return Standard_False;
1736   }
1737
1738   aProjector.LowerDistanceParameters(ULD, VLD);
1739   aP2D.SetCoord(ULD, VLD);
1740   //
1741   bFlag=aContext->IsPointInOnFace (aF, aP2D);
1742   //
1743   if (!bFlag) {
1744     return bFlag;
1745   }
1746   return bFlag;
1747 }
1748
1749 //=======================================================================
1750 //function : IsMicroEdge
1751 //purpose  : 
1752 //=======================================================================
1753   Standard_Boolean BOPTools_AlgoTools::IsMicroEdge(const TopoDS_Edge& aE,
1754                                                    const Handle(BOPInt_Context)& aCtx) 
1755 {
1756   Standard_Boolean bRet;
1757   Standard_Integer iErr;
1758   Standard_Real aT1, aT2, aTmp;
1759   Handle(Geom_Curve) aC3D;
1760   TopoDS_Vertex aV1, aV2;
1761   //
1762   bRet=(BRep_Tool::Degenerated(aE) ||
1763         !BRep_Tool::IsGeometric(aE));
1764   if (bRet) {
1765     return bRet;
1766   }
1767   //
1768   aC3D=BRep_Tool::Curve(aE, aT1, aT2);
1769   TopExp::Vertices(aE, aV1, aV2);
1770   aT1=BRep_Tool::Parameter(aV1, aE);
1771   aT2=BRep_Tool::Parameter(aV2, aE);
1772   if (aT2<aT1) {
1773     aTmp=aT1;
1774     aT1=aT2;
1775     aT2=aTmp;
1776   }
1777   //
1778   BOPInt_ShrunkRange aSR;
1779   aSR.SetData(aE, aT1, aT2, aV1, aV2, aCtx);
1780   aSR.Perform();
1781   iErr=aSR.ErrorStatus();
1782   bRet = !(iErr==0);
1783   //
1784   return bRet;
1785 }
1786
1787 //=======================================================================
1788 //function : GetFaceDir
1789 //purpose  : Get binormal direction for the face in the point aP
1790 //=======================================================================
1791   void GetFaceDir(const TopoDS_Edge& aE,
1792                   const TopoDS_Face& aF,
1793                   const gp_Pnt& aP,
1794                   const Standard_Real aT,
1795                   const gp_Dir& aDTgt,
1796                   gp_Dir& aDN,
1797                   gp_Dir& aDB,
1798                   Handle(BOPInt_Context)& theContext)
1799 {
1800   BOPTools_AlgoTools3D::GetNormalToFaceOnEdge(aE, aF, aT, aDN);
1801   if (aF.Orientation()==TopAbs_REVERSED){
1802     aDN.Reverse();
1803   }
1804   aDB = aDN^aDTgt;
1805   //
1806   Handle(Geom_Surface) aS = BRep_Tool::Surface(aF);
1807   GeomAdaptor_Surface aGAS(aS);
1808   if (aGAS.GetType()!=GeomAbs_Plane) {
1809     gp_Pnt aPx;
1810     if (!FindPointInFace(aE, aF, aP, aT, aDB, aPx, theContext, aGAS)) {
1811       gp_Pnt2d aPx2D;
1812       Handle(Geom_Plane) aPL;
1813       GeomAPI_ProjectPointOnSurf aProj;
1814       //
1815       BOPTools_AlgoTools3D::PointNearEdge(aE, aF, aT, aPx2D, aPx, theContext);
1816       aPL = new Geom_Plane(aP, aDTgt);
1817       aProj.Init(aPx, aPL);
1818       aPx = aProj.NearestPoint();
1819       gp_Vec aVec(aP, aPx);
1820       aDB.SetXYZ(aVec.XYZ());
1821     }
1822   }
1823 }
1824
1825 //=======================================================================
1826 //function : FindPointInFace
1827 //purpose  : Find a point in the face in direction of <aDB>
1828 //=======================================================================
1829   Standard_Boolean FindPointInFace(const TopoDS_Edge& aE,
1830                                    const TopoDS_Face& aF,
1831                                    const gp_Pnt& aP,
1832                                    const Standard_Real /*aT*/,
1833                                    gp_Dir& aDB,
1834                                    gp_Pnt& aPOut,
1835                                    Handle(BOPInt_Context)& theContext,
1836                                    const GeomAdaptor_Surface& aGAS) 
1837 {
1838   Standard_Integer aNbItMax;
1839   Standard_Real aDt, aDtMin, aTolE, aTolF, aDist;
1840   Standard_Boolean bRet;
1841   gp_Pnt aP1;
1842   //
1843   bRet = Standard_False;
1844   aTolE = BRep_Tool::Tolerance(aE);
1845   aTolF = BRep_Tool::Tolerance(aF);
1846   aDt = 2*(aTolE+aTolF);
1847   //
1848   aDtMin=5.e-4;
1849   if (aDt < aDtMin) {
1850     Standard_Real aR;
1851     GeomAbs_SurfaceType aSType=aGAS.GetType();
1852     switch (aSType) {
1853     case GeomAbs_Cylinder:
1854       aR = aGAS.Cylinder().Radius();
1855       break;
1856     case GeomAbs_Cone:
1857       aR = aGAS.Cone().RefRadius();
1858       break;
1859     case GeomAbs_Sphere:
1860       aR = aGAS.Sphere().Radius();
1861       break;
1862     case GeomAbs_Torus:
1863       aR = aGAS.Torus().MinorRadius();
1864       break;
1865     case GeomAbs_SurfaceOfRevolution:
1866       aR=1.;
1867       break;
1868     default:
1869       aR=0.001;
1870     }
1871     if (aR < 0.01) {
1872       aDtMin=5.e-6;
1873     }
1874     else if (aR > 100.) {
1875       aDtMin*=10;
1876     }
1877     if (aDt < aDtMin) {
1878       aDt=aDtMin;
1879     }
1880   }
1881   //
1882   GeomAPI_ProjectPointOnSurf& aProj=theContext->ProjPS(aF);
1883   aNbItMax = 15;
1884   //
1885   do {
1886     aP1.SetCoord(aP.X()+aDt*aDB.X(),
1887                  aP.Y()+aDt*aDB.Y(),
1888                  aP.Z()+aDt*aDB.Z());
1889     //
1890     aProj.Perform(aP1);
1891     if (!aProj.IsDone()) {
1892       return bRet;
1893     }
1894     aPOut = aProj.NearestPoint();
1895     aDist = aProj.LowerDistance();
1896     //
1897     gp_Vec aV(aP, aPOut);
1898     aDB.SetXYZ(aV.XYZ());
1899   } while (aDist>Precision::Angular() && --aNbItMax);
1900   //
1901   bRet = aDist < Precision::Angular();
1902   return bRet;
1903 }