0023676: Intersection between edge and a face produces extra vertex
[occt.git] / src / IntTools / IntTools_EdgeFace.cxx
1 // Created on: 2001-02-26
2 // Created by: Peter KURNEV
3 // Copyright (c) 2001-2012 OPEN CASCADE SAS
4 //
5 // The content of this file is subject to the Open CASCADE Technology Public
6 // License Version 6.5 (the "License"). You may not use the content of this file
7 // except in compliance with the License. Please obtain a copy of the License
8 // at http://www.opencascade.org and read it completely before using this file.
9 //
10 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
11 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
12 //
13 // The Original Code and all software distributed under the License is
14 // distributed on an "AS IS" basis, without warranty of any kind, and the
15 // Initial Developer hereby disclaims all such warranties, including without
16 // limitation, any warranties of merchantability, fitness for a particular
17 // purpose or non-infringement. Please see the License for the specific terms
18 // and conditions governing the rights and limitations under the License.
19
20
21 #include <IntTools_EdgeFace.ixx>
22
23
24
25 #include <IntTools_CArray1OfReal.hxx>
26 #include <IntTools.hxx>
27 #include <IntTools_CArray1OfInteger.hxx>
28 #include <IntTools_Range.hxx>
29 #include <IntTools_Tools.hxx>
30 #include <IntTools_Array1OfRange.hxx>
31 #include <IntTools_QuickSortRange.hxx>
32 #include <IntTools_CompareRange.hxx>
33 #include <IntTools_CommonPrt.hxx>
34 #include <IntTools_Root.hxx>
35 #include <IntTools_BeanFaceIntersector.hxx>
36 #include <IntTools_Context.hxx>
37
38 #include <BRep_Tool.hxx>
39
40 #include <GeomAdaptor_Surface.hxx>
41 #include <GeomAdaptor_Curve.hxx>
42
43 #include <Geom_Surface.hxx>
44 #include <Geom_Curve.hxx>
45
46 #include <GeomAPI_ProjectPointOnSurf.hxx>
47
48 #include <Precision.hxx>
49
50 #include <Bnd_Box.hxx>
51 #include <BndLib_AddSurface.hxx>
52
53 #include <gp_Cylinder.hxx>
54 #include <gp_Ax1.hxx>
55 #include <gp_Lin.hxx>
56 #include <gp_Cone.hxx>
57 #include <gp_Torus.hxx>
58 #include <gp_Circ.hxx>
59 #include <gp_Pln.hxx>
60
61
62 #include <Extrema_ExtCS.hxx>
63 #include <Extrema_POnCurv.hxx>
64 #include <Extrema_POnSurf.hxx>
65
66 // modified by NIZHNY-MKK  Thu Jul 21 11:35:59 2005
67 #include <IntCurveSurface_HInter.hxx>
68 #include <GeomAdaptor_HCurve.hxx>
69 #include <GeomAdaptor_HSurface.hxx>
70 #include <IntCurveSurface_IntersectionPoint.hxx>
71
72 #ifdef WNT
73 #pragma warning ( disable : 4101 )
74 #endif
75
76 static
77   Standard_Boolean IsCoplanar (const BRepAdaptor_Curve&  ,
78                                const BRepAdaptor_Surface& );
79 static
80   Standard_Boolean IsRadius (const BRepAdaptor_Curve& aCurve ,
81                              const BRepAdaptor_Surface& aSurface);
82 static
83   Standard_Integer AdaptiveDiscret (const Standard_Integer iDiscret,
84                                     const BRepAdaptor_Curve& aCurve ,
85                                     const BRepAdaptor_Surface& aSurface);
86
87 //=======================================================================
88 //function : IntTools_EdgeFace::IntTools_EdgeFace
89 //purpose  : 
90 //=======================================================================
91   IntTools_EdgeFace::IntTools_EdgeFace()
92 {
93   myTolE=1.e-7;
94   myTolF=1.e-7;
95   myDiscret=30;
96   myEpsT   =1e-12;
97   myEpsNull=1e-12;
98   myDeflection=0.01;
99   myIsDone=Standard_False;
100   myErrorStatus=1;
101   myParallel=Standard_False;
102   myPar1=0.;
103 }
104 //=======================================================================
105 //function : SetContext
106 //purpose  : 
107 //=======================================================================
108 void IntTools_EdgeFace::SetContext(const Handle(IntTools_Context)& theContext) 
109 {
110   myContext = theContext;
111 }
112
113 //=======================================================================
114 //function : Context
115 //purpose  : 
116 //=======================================================================
117 const Handle(IntTools_Context)& IntTools_EdgeFace::Context()const 
118 {
119   return myContext;
120 }
121 //=======================================================================
122 //function : SetEdge
123 //purpose  : 
124 //=======================================================================
125   void IntTools_EdgeFace::SetEdge(const TopoDS_Edge& anEdge)
126 {
127   myEdge=anEdge;
128 }
129
130 //=======================================================================
131 //function : SetFace
132 //purpose  : 
133 //=======================================================================
134   void IntTools_EdgeFace::SetFace(const TopoDS_Face& aFace)
135 {
136   myFace=aFace;
137 }
138
139 //=======================================================================
140 //function : SetTolE
141 //purpose  : 
142 //=======================================================================
143   void IntTools_EdgeFace::SetTolE(const Standard_Real aTol) 
144 {
145   myTolE=aTol;
146
147 //=======================================================================
148 //function : SetTolF
149 //purpose  : 
150 //=======================================================================
151   void IntTools_EdgeFace::SetTolF(const Standard_Real aTol) 
152 {
153   myTolF=aTol;
154
155
156 //=======================================================================
157 //function : SetDiscretize
158 //purpose  : 
159 //=======================================================================
160   void IntTools_EdgeFace::SetDiscretize(const Standard_Integer aDiscret)
161 {
162   myDiscret=aDiscret;
163 }
164 //=======================================================================
165 //function : SetDeflection
166 //purpose  : 
167 //=======================================================================
168   void IntTools_EdgeFace::SetDeflection(const Standard_Real aDefl) 
169 {
170   myDeflection=aDefl;
171
172 //=======================================================================
173 //function : SetEpsilonT
174 //purpose  : 
175 //=======================================================================
176   void IntTools_EdgeFace::SetEpsilonT(const Standard_Real anEpsT) 
177 {
178   myEpsT=anEpsT;
179
180 //=======================================================================
181 //function : SetEpsilonNull
182 //purpose  : 
183 //=======================================================================
184   void IntTools_EdgeFace::SetEpsilonNull(const Standard_Real anEpsNull) 
185 {
186   myEpsNull=anEpsNull;
187
188
189 //=======================================================================
190 //function : SetRange
191 //purpose  : 
192 //=======================================================================
193   void IntTools_EdgeFace::SetRange(const Standard_Real aFirst,
194                                    const Standard_Real aLast) 
195 {
196   myRange.SetFirst (aFirst);
197   myRange.SetLast  (aLast);
198
199
200 //=======================================================================
201 //function : SetRange
202 //purpose  : 
203 //=======================================================================
204   void IntTools_EdgeFace::SetRange(const IntTools_Range& aRange) 
205 {
206   myRange.SetFirst (aRange.First());
207   myRange.SetLast  (aRange.Last());
208
209 //=======================================================================
210 //function : IsDone
211 //purpose  : 
212 //=======================================================================
213 Standard_Boolean IntTools_EdgeFace::IsDone()const 
214 {
215   return myIsDone;
216
217 //=======================================================================
218 //function : ErrorStatus
219 //purpose  : 
220 //=======================================================================
221 Standard_Integer IntTools_EdgeFace::ErrorStatus()const 
222 {
223   return myErrorStatus;
224 }
225 //=======================================================================
226 //function : CommonParts
227 //purpose  : 
228 //=======================================================================
229 const IntTools_SequenceOfCommonPrts& IntTools_EdgeFace::CommonParts() const 
230 {
231   return mySeqOfCommonPrts;
232 }
233 //=======================================================================
234 //function : Range
235 //purpose  : 
236 //=======================================================================
237 const IntTools_Range&  IntTools_EdgeFace::Range() const
238 {
239   return myRange;
240
241
242 //=======================================================================
243 //function : CheckData
244 //purpose  : 
245 //=======================================================================
246 void IntTools_EdgeFace::CheckData()
247 {
248   if (BRep_Tool::Degenerated(myEdge)) {
249     myErrorStatus=2;
250   }
251   if (!BRep_Tool::IsGeometric(myEdge)) { 
252      myErrorStatus=3;
253   }
254 }
255 //=======================================================================
256 //function : Prepare
257 //purpose  : 
258 //=======================================================================
259   void IntTools_EdgeFace::Prepare() 
260 {
261   Standard_Integer pri;
262   IntTools_CArray1OfReal aPars;
263  
264   //
265   // 1.Prepare Curve's data and Surface's data
266   myC.Initialize(myEdge);
267   GeomAbs_CurveType aCurveType;
268   aCurveType=myC.GetType();
269   //
270   // 2.Prepare myCriteria
271   if (aCurveType==GeomAbs_BSplineCurve||
272         aCurveType==GeomAbs_BezierCurve) {
273     myCriteria=1.5*myTolE+myTolF;
274   }
275   else {
276     myCriteria=myTolE+myTolF;
277   }
278   // 2.a myTmin, myTmax
279   myTmin=myRange.First();
280   myTmax=myRange.Last();
281   // 2.b myFClass2d
282   myS.Initialize (myFace,Standard_True);
283   myFClass2d.Init(myFace, 1.e-6);
284   //
285   // 2.c Prepare adaptive myDiscret
286   myDiscret=AdaptiveDiscret(myDiscret, myC, myS);
287   //
288   //
289   // 3.Prepare myPars 
290   pri = IntTools::PrepareArgs(myC, myTmax, myTmin, myDiscret, myDeflection, aPars);
291   if (pri) {
292     myErrorStatus=6;
293     return;
294   }
295   // 4.
296   //ProjectableRanges
297   Standard_Integer i, iProj, aNb, aNbProj, ind0, ind1;
298   Standard_Real t0, t1, tRoot;
299   
300   //
301   // Table of Projection's function values
302   aNb=aPars.Length();
303   IntTools_CArray1OfInteger anArrProjectability;
304   anArrProjectability.Resize(aNb);
305   
306   for (iProj=0, i=0; i<aNb; i++) {
307     t0=aPars(i);
308     aNbProj=IsProjectable (t0); 
309     
310     anArrProjectability(i)=0;
311     if (aNbProj) {
312       anArrProjectability(i)=1;
313       iProj++;
314     }
315   }
316   //
317   // Checking
318   if (!iProj ) {
319     myErrorStatus=7;
320     return;
321   }
322   
323   //
324   // Projectable Ranges
325   IntTools_Range aRange;
326   
327   ind0=anArrProjectability(0);
328   if (ind0) {
329     t0=aPars(0);
330     aRange.SetFirst(t0);
331   }
332   
333   for(i=1; i<aNb; i++) {
334     ind1=anArrProjectability(i);
335     t0=aPars(i-1);
336     t1=aPars(i);
337
338     if (i==(aNb-1)) {
339       if (ind1 && ind0) {
340         aRange.SetLast(t1);
341         myProjectableRanges.Append(aRange);
342       }
343       if (ind1 && !ind0) {
344         FindProjectableRoot(t0, t1, ind0, ind1, tRoot);
345         aRange.SetFirst(tRoot);
346         aRange.SetLast(t1);
347         myProjectableRanges.Append(aRange);
348       }
349       //
350       if (ind0 && !ind1) {
351         FindProjectableRoot(t0, t1, ind0, ind1, tRoot);
352         aRange.SetLast(tRoot);
353         myProjectableRanges.Append(aRange);
354       }
355       //
356       break;
357     }
358     
359     if (ind0 != ind1) {
360       FindProjectableRoot(t0, t1, ind0, ind1, tRoot);
361       
362       if (ind0 && !ind1) {
363         aRange.SetLast(tRoot);
364         myProjectableRanges.Append(aRange);
365       }
366       else {
367         aRange.SetFirst(tRoot);
368       }
369     } // if (ind0 != ind1)
370     ind0=ind1;
371   } // for(i=1; i<aNb; i++) {
372 }
373
374 //=======================================================================
375 //function : FindProjectableRoot
376 //purpose  : 
377 //=======================================================================
378   void IntTools_EdgeFace::FindProjectableRoot (const Standard_Real tt1,
379                                                const Standard_Real tt2,
380                                                const Standard_Integer ff1,
381                                                const Standard_Integer ff2,
382                                                Standard_Real& tRoot)
383 {
384   Standard_Real tm, t1, t2, aEpsT;
385   Standard_Integer anIsProj1, anIsProj2, anIsProjm;
386   aEpsT=0.5*myEpsT;
387   //
388   // Root is inside [tt1, tt2]
389   t1=tt1;  
390   t2=tt2;
391   anIsProj1=ff1;
392   anIsProj2=ff2;
393   
394   while (1) {
395     if (fabs(t1-t2) < aEpsT) {
396       tRoot=(anIsProj1) ? t1 : t2;
397       return;
398     }
399     tm=.5*(t1+t2);
400     anIsProjm=IsProjectable(tm);
401     
402     if (anIsProjm != anIsProj1) {
403       t2=tm;
404       anIsProj2=anIsProjm;
405     }
406     else {
407       t1=tm;
408       anIsProj1=anIsProjm;
409     }
410   }
411 }
412 //=======================================================================
413 //function : IsProjectable
414 //purpose  : 
415 //=======================================================================
416   Standard_Boolean IntTools_EdgeFace::IsProjectable(const Standard_Real t) const
417 {
418   Standard_Boolean bFlag;
419   Standard_Real Umin, Usup, Vmin, Vsup;
420
421   Umin=myS.FirstUParameter();
422   Usup=myS.LastUParameter();
423   Vmin=myS.FirstVParameter();
424   Vsup=myS.LastVParameter ();
425   
426   gp_Pnt P;
427   myC.D0(t, P);
428   GeomAPI_ProjectPointOnSurf aProjector;
429   //
430   Standard_Real ULD, VLD;
431
432   GeomAPI_ProjectPointOnSurf& aLocProj = myContext->ProjPS(myFace);
433   aLocProj.Perform(P);
434   bFlag = aLocProj.IsDone();
435   
436   if(bFlag) {
437     aLocProj.LowerDistanceParameters(ULD, VLD);
438   }
439   //
440
441   if (bFlag) {
442     bFlag=Standard_False;
443
444     // 
445     TopAbs_State aState;
446     gp_Pnt2d aP2d(ULD, VLD);
447
448     aState = myContext->FClass2d(myFace).Perform(aP2d);
449     //
450     
451     if (aState==TopAbs_IN || aState==TopAbs_ON) {
452       bFlag=Standard_True;
453     }
454   }
455   return bFlag;
456 }
457
458 //=======================================================================
459 //function : DistanceFunction
460 //purpose  : 
461 //=======================================================================
462   Standard_Real IntTools_EdgeFace::DistanceFunction(const Standard_Real t)
463 {
464   Standard_Real Umin, Usup, Vmin, Vsup, aD;
465   //
466   gp_Pnt P;
467   myC.D0(t, P);
468   //
469   Standard_Boolean bIsEqDistance;
470
471   bIsEqDistance= IntTools_EdgeFace::IsEqDistance(P, myS, 1.e-7, aD); 
472   if (bIsEqDistance) {
473     aD=aD-myCriteria;
474     return aD; 
475   }
476   
477   Umin=myS.FirstUParameter();
478   Usup=myS.LastUParameter();
479   Vmin=myS.FirstVParameter();
480   Vsup=myS.LastVParameter ();
481   
482   //
483   Standard_Boolean bFlag = Standard_False;
484
485   GeomAPI_ProjectPointOnSurf& aLocProj = myContext->ProjPS(myFace);
486   aLocProj.Perform(P);
487   bFlag = aLocProj.IsDone();
488   
489   if(bFlag) {
490     aD = aLocProj.LowerDistance();
491   }
492   //
493   
494   if (!bFlag) {
495     myErrorStatus=11;
496     return 99.;
497   }
498   
499   // 
500   //   aD=aProjector.LowerDistance();
501   // 
502   aD=aD-myCriteria;
503   return aD; 
504 }
505 //
506 //=======================================================================
507 //function : IsEqDistance
508 //purpose  : 
509 //=======================================================================
510   Standard_Boolean IntTools_EdgeFace::IsEqDistance(const gp_Pnt& aP,
511                                                    const BRepAdaptor_Surface& aBAS,
512                                                    const Standard_Real aTol,
513                                                    Standard_Real& aD)
514 {
515   Standard_Boolean bRetFlag=Standard_True;
516
517   GeomAbs_SurfaceType aSurfType=aBAS.GetType();
518
519   if (aSurfType==GeomAbs_Cylinder) {
520     gp_Cylinder aCyl=aBAS.Cylinder();
521     const gp_Ax1& anAx1  =aCyl.Axis();
522     gp_Lin aLinAxis(anAx1);
523     Standard_Real aDC, aRadius=aCyl.Radius();
524     aDC=aLinAxis.Distance(aP);
525     if (aDC < aTol) {
526       aD=aRadius;
527       return bRetFlag; 
528     }
529   }
530
531   if (aSurfType==GeomAbs_Cone) {
532     gp_Cone aCone=aBAS.Cone();
533     const gp_Ax1& anAx1  =aCone.Axis();
534     gp_Lin aLinAxis(anAx1);
535     Standard_Real aDC, aRadius, aDS, aSemiAngle;
536     aDC=aLinAxis.Distance(aP);
537     if (aDC < aTol) {
538       gp_Pnt anApex=aCone.Apex();
539       aSemiAngle=aCone.SemiAngle();
540       aDS=aP.Distance(anApex);
541       
542       aRadius=aDS*tan(aSemiAngle);
543       aD=aRadius;
544       return bRetFlag; 
545     }
546   }
547
548   if (aSurfType==GeomAbs_Torus) {
549     Standard_Real aMajorRadius, aMinorRadius, aDC;
550
551     gp_Torus aTorus=aBAS.Torus();
552     gp_Pnt aPLoc=aTorus.Location();
553     aMajorRadius=aTorus.MajorRadius();
554     
555     aDC=fabs(aPLoc.Distance(aP)-aMajorRadius);
556     if (aDC < aTol) {
557       aMinorRadius=aTorus.MinorRadius();
558       aD=aMinorRadius;
559       return bRetFlag; 
560     }
561   }
562   return !bRetFlag; 
563 }
564 //
565 //=======================================================================
566 //function : PrepareArgsFuncArrays
567 //purpose  : Obtain 
568 //           myFuncArray and myArgsArray for the interval [ta, tb]
569 //=======================================================================               
570   void IntTools_EdgeFace::PrepareArgsFuncArrays(const Standard_Real ta,
571                                                 const Standard_Real tb)
572 {
573   IntTools_CArray1OfReal anArgs, aFunc;
574   Standard_Integer i, aNb, pri;
575   Standard_Real t, f, f1;
576   //
577   // Prepare values of arguments for the interval [ta, tb]
578   pri=IntTools::PrepareArgs (myC, tb, ta, myDiscret, myDeflection, anArgs);
579   
580   if (pri) {
581     myErrorStatus=8;
582     return;
583   }
584   //...
585   aNb=anArgs.Length();
586
587   if (!aNb){
588     myErrorStatus=9;
589     return;
590   }
591   //
592   // Prepare values of functions for the interval [ta, tb]
593   aFunc.Resize(aNb);
594   for (i=0; i<aNb; i++) {
595     t=anArgs(i);
596     f1=DistanceFunction(t);
597     f=f1+myCriteria; 
598
599     if (myErrorStatus==11)
600       return;
601     
602     if (f1 < myEpsNull) { 
603       f=0.;
604     }
605     aFunc(i)=f;
606   }
607   //
608   // Add points where the derivative = 0  
609   AddDerivativePoints(anArgs, aFunc);
610
611 }
612
613
614 //=======================================================================
615 //function : AddDerivativePoints
616 //purpose  : 
617 //=======================================================================
618   void IntTools_EdgeFace::AddDerivativePoints(const IntTools_CArray1OfReal& t,
619                                               const IntTools_CArray1OfReal& f)  
620 {
621   Standard_Integer i, j, n, k, nn=100;
622   Standard_Real fr, tr, tr1, dEpsNull=10.*myEpsNull;
623   IntTools_CArray1OfReal fd;
624   TColStd_SequenceOfReal aTSeq, aFSeq;  
625
626   n=t.Length();
627   fd.Resize(n+1);
628   //
629   // Table of derivatives
630   Standard_Real dfx, tx, tx1, fx, fx1, dt=1.e-6;
631   // Left limit
632   tx=t(0);
633   tx1=tx+dt;
634   fx=f(0);
635   fx1=DistanceFunction(tx1);
636   fx1=fx1+myCriteria;
637   if (fx1 < myEpsNull) { 
638     fx1=0.;
639   }
640   dfx=(fx1-fx)/dt;
641   fd(0)=dfx;
642   
643   if (fabs(fd(0)) < dEpsNull){
644     fd(0)=0.;
645   }
646   
647
648   k=n-1;
649   for (i=1; i<k; i++) {
650     Standard_Real ti, ti1;
651     ti=t(i);
652     ti1=t(i-1);
653     fd(i)=.5*(f(i+1)-f(i-1))/(t(i)-t(i-1));
654     if (fabs(fd(i)) < dEpsNull){
655       fd(i)=0.;
656     }
657   }
658   // Right limit
659   tx=t(n-1);
660   tx1=tx-dt;
661   fx=f(n-1);
662   fx1=DistanceFunction(tx1);
663   fx1=fx1+myCriteria;
664   if (fx1 < myEpsNull) { 
665     fx1=0.;
666   }
667   dfx=(fx-fx1)/dt;
668   fd(n-1)=dfx;
669   
670   if (fabs(fd(n-1)) < dEpsNull){
671     fd(n-1)=0.;
672   }
673   //
674   // Finding the range where the derivatives have different signs
675   // for neighbouring points
676   for (i=1; i<n; i++) {
677     Standard_Real fd1, fd2, t1, t2;
678     t1 =t(i-1);
679     t2 =t(i);
680     fd1=fd(i-1);
681     fd2=fd(i);
682
683     if (fd1*fd2 < 0.) {
684       if (fabs(fd1) < myEpsNull) {
685         tr=t1;
686         fr=DistanceFunction(tr);//fd1;
687       }
688       else if (fabs(fd2) < myEpsNull) {
689         tr=t2;
690         fr=DistanceFunction(tr);
691       }
692       else {
693         tr=FindSimpleRoot(2, t1, t2, fd1);
694         fr=DistanceFunction(tr);
695       }
696       
697       aTSeq.Append(tr);
698       aFSeq.Append(fr);
699     }
700   } // end of for (i=1; i<n; i++)
701   //
702   // remove identical t, f
703   nn=aTSeq.Length();
704   if (nn) {
705     for (i=1; i<=aTSeq.Length(); i++) {
706       tr=aTSeq(i);
707       for (j=0; j<n; j++) {
708         tr1=t(j);
709         if (fabs (tr1-tr) < myEpsT) {
710           aTSeq.Remove(i);
711           aFSeq.Remove(i);
712         }
713       }
714     }
715     nn=aTSeq.Length();
716   }
717   //
718   // sorting args and funcs in increasing order
719   if (nn) {
720     k=nn+n;
721     IntTools_Array1OfRange anArray1OfRange(1, k);
722     for (i=1; i<=n; i++) {
723       anArray1OfRange(i).SetFirst(t(i-1));
724       anArray1OfRange(i).SetLast (f(i-1));
725     }
726     for (i=1; i<=nn; i++) {
727       anArray1OfRange(n+i).SetFirst(aTSeq(i));
728       anArray1OfRange(n+i).SetLast (aFSeq(i));
729     }
730     
731     IntTools_QuickSortRange aQuickSortRange;
732     IntTools_CompareRange aComparator;
733     aQuickSortRange.Sort (anArray1OfRange, aComparator);
734     
735     // filling the  output arrays
736     myArgsArray.Resize(k);
737     myFuncArray.Resize(k);
738     for (i=1; i<=k; i++) {
739       myArgsArray(i-1)=anArray1OfRange(i).First();
740       myFuncArray(i-1)=anArray1OfRange(i).Last ();
741     }
742   }
743   
744   else { // nn=0
745     myArgsArray.Resize(n);
746     myFuncArray.Resize(n);
747     for (i=0; i<n; i++) {
748       myArgsArray(i)=t(i); 
749       myFuncArray(i)=f(i); 
750     }
751   }
752 }
753
754 //=======================================================================
755 //function : DerivativeFunction
756 //purpose  : 
757 //=======================================================================
758   Standard_Real IntTools_EdgeFace::DerivativeFunction(const Standard_Real t2)
759 {
760   Standard_Real t1, t3, aD1, aD2, aD3;
761   Standard_Real dt=1.e-9;
762   t1=t2-dt;
763   aD1=DistanceFunction(t1);
764   t3=t2+dt;
765   aD3=DistanceFunction(t3);
766   
767   aD2=.5*(aD3-aD1)/dt;
768   return aD2; 
769 }
770
771 //=======================================================================
772 //function : FindSimpleRoot
773 //purpose  : [private]
774 //=======================================================================
775   Standard_Real IntTools_EdgeFace::FindSimpleRoot (const Standard_Integer IP,
776                                                    const Standard_Real tA,
777                                                    const Standard_Real tB,
778                                                    const Standard_Real fA)
779 {
780   Standard_Real r, a, b, y, x0, s;
781   
782   a=tA; b=tB; r=fA;
783   
784   while (1) {
785     x0=.5*(a+b);
786
787     if (IP==1)
788       y=DistanceFunction(x0);
789     else 
790       y=DerivativeFunction(x0);
791     
792     if (fabs(b-a) < myEpsT || y==0.) {
793       return x0;
794     }
795     
796         
797     s=y*r;
798
799     if (s<0.) {
800       b=x0;
801       continue;
802     }
803
804     if (s>0.) {
805       a=x0; r=y;
806     }
807   }
808 }
809 //=======================================================================
810 //function : FindGoldRoot
811 //purpose  : [private]
812 //=======================================================================
813   Standard_Real IntTools_EdgeFace::FindGoldRoot (const Standard_Real tA,
814                                                  const Standard_Real tB,
815                                                  const Standard_Real coeff)
816 {
817   Standard_Real gs=0.61803399;
818   Standard_Real a, b, xp, xl, yp, yl;
819
820   a=tA;  b=tB;
821   
822   xp=a+(b-a)*gs;
823   xl=b-(b-a)*gs;
824   yp=coeff*DistanceFunction(xp);
825   yl=coeff*DistanceFunction(xl);
826   
827  
828   while (1) {
829     
830     if (fabs(b-a) < myEpsT) {
831       return .5*(b+a);
832     }
833     
834     if (yp < yl) {
835       a=xl;
836       xl=xp;
837       xp=a+(b-a)*gs;
838       yp=coeff*DistanceFunction(xp);
839     }
840     
841     else {
842       b=xp;
843       xp=xl;
844       yp=yl;
845       xl=b-(b-a)*gs;
846       yl=coeff*DistanceFunction(xl);
847     }
848   }
849 }  
850
851 //=======================================================================
852 //function : MakeType
853 //purpose  : 
854 //=======================================================================
855   Standard_Integer IntTools_EdgeFace::MakeType(IntTools_CommonPrt&  aCommonPrt)
856 {
857   Standard_Real  af1, al1;
858   Standard_Real dt, df1, df2, tm;
859   Standard_Boolean bAllNullFlag;
860   //
861   bAllNullFlag=aCommonPrt.AllNullFlag();
862   if (bAllNullFlag) {
863     aCommonPrt.SetType(TopAbs_EDGE);
864     return 0;
865   }
866   //
867   aCommonPrt.Range1(af1, al1);
868
869   {
870     gp_Pnt aPF, aPL;
871     myC.D0(af1, aPF);
872     myC.D0(al1, aPL);
873     df1=aPF.Distance(aPL);
874     Standard_Boolean isWholeRange = Standard_False;
875     
876     if((Abs(af1 - myRange.First()) < myC.Resolution(myCriteria)) &&
877        (Abs(al1 - myRange.Last()) < myC.Resolution(myCriteria)))
878       isWholeRange = Standard_True;
879
880     
881     if ((df1 > myCriteria * 2.) && isWholeRange) {
882       aCommonPrt.SetType(TopAbs_EDGE);
883     }
884     else {
885       if(isWholeRange) {
886         tm = (af1 + al1) * 0.5;
887
888         if(aPF.Distance(myC.Value(tm)) > myCriteria * 2.) {
889           aCommonPrt.SetType(TopAbs_EDGE);
890           return 0;
891         }
892       }
893
894       if(!CheckTouch(aCommonPrt, tm)) {
895         tm = (af1 + al1) * 0.5;
896       }
897       aCommonPrt.SetType(TopAbs_VERTEX);
898       aCommonPrt.SetVertexParameter1(tm);
899       aCommonPrt.SetRange1 (af1, al1);
900     }
901          return 0;
902   }
903   //
904   /*
905   dt=al1-af1;
906   if (dt<1.e-5) {
907     gp_Pnt aPF, aPL;
908     myC.D0(af1, aPF);
909     myC.D0(al1, aPL);
910     df1=aPF.Distance(aPL);
911     if (df1<myCriteria) {
912       //
913       tm=.5*(af1+al1);
914       aCommonPrt.SetType(TopAbs_VERTEX);
915       aCommonPrt.SetVertexParameter1(tm);
916       aCommonPrt.SetRange1 (af1, al1);
917       return 0;
918     }
919   }
920   //
921   IsIntersection (af1, al1);
922   //
923   if (!myParallel) {
924     aCommonPrt.SetType(TopAbs_VERTEX);
925     aCommonPrt.SetVertexParameter1(myPar1);
926     aCommonPrt.SetRange1 (af1, al1);
927   }
928   else {
929     dt=al1-af1;
930     if (dt<1.e-5) {
931       df1=DistanceFunction(af1);
932       df2=DistanceFunction(al1);
933       tm=(df1 < df2) ? af1 : al1;
934       aCommonPrt.SetType(TopAbs_VERTEX);
935       aCommonPrt.SetVertexParameter1(tm);
936       aCommonPrt.SetRange1 (af1, al1);
937     }
938
939     else {
940       aCommonPrt.SetType(TopAbs_EDGE);
941     }
942   }
943   */
944   return 0;
945 }
946
947
948 //=======================================================================
949 //function : IsIntersection
950 //purpose  : 
951 //=======================================================================
952   void IntTools_EdgeFace::IsIntersection (const Standard_Real ta, 
953                                           const Standard_Real tb) 
954 {
955   IntTools_CArray1OfReal anArgs, aFunc;
956   Standard_Integer i, aNb, pri, aCnt=0;
957   //
958   Standard_Integer aCntIncreasing=1, aCntDecreasing=1;
959   Standard_Real t, f, f1;
960   //
961   // Prepare values of arguments for the interval [ta, tb]
962   pri=IntTools::PrepareArgs (myC, tb, ta, myDiscret, myDeflection, anArgs);
963   aNb=anArgs.Length();
964   
965   aFunc.Resize(aNb);
966   for (i=0; i<aNb; i++) {
967     t=anArgs(i);
968     
969     f1=DistanceFunction(t);
970     f=f1+myCriteria; 
971
972     if (fabs(f1) < myEpsNull) { 
973       aCnt++;
974       f=0.;
975     }
976     aFunc(i)=f;
977     //
978     if (i) {
979       if (aFunc(i)>aFunc(i-1)) {
980         aCntIncreasing++;
981       }
982       if (aFunc(i)<aFunc(i-1)) {
983         aCntDecreasing++;
984       }
985     }
986     //
987   }
988
989   if (aCnt==aNb) {
990     myParallel=Standard_True;
991     return;
992   }
993   
994   FindDerivativeRoot(anArgs, aFunc);
995
996   //
997   if (myParallel) {
998     if (!(myC.GetType()==GeomAbs_Line 
999           && 
1000           myS.GetType()==GeomAbs_Cylinder)) {
1001       if (aCntDecreasing==aNb) {
1002         myPar1=anArgs(aNb-1);
1003         myParallel=Standard_False;
1004       }
1005       if (aCntIncreasing==aNb) {
1006         myPar1=anArgs(0);
1007         myParallel=Standard_False;
1008       }
1009     }
1010   }
1011   //
1012   return ;
1013 }
1014
1015 //=======================================================================
1016 //function : FindDerivativeRoot
1017 //purpose  : 
1018 //=======================================================================
1019   void IntTools_EdgeFace::FindDerivativeRoot(const IntTools_CArray1OfReal& t,
1020                                              const IntTools_CArray1OfReal& f)  
1021 {
1022   Standard_Integer i, n, k;
1023   Standard_Real fr, tr;
1024   IntTools_CArray1OfReal fd;
1025   TColStd_SequenceOfReal aTSeq, aFSeq;  
1026   
1027   myPar1=0.;
1028   myParallel=Standard_True;
1029   
1030   n=t.Length();
1031   fd.Resize(n+1);
1032   //
1033   // Table of derivatives
1034   fd(0)=(f(1)-f(0))/(t(1)-t(0));
1035   if (fabs(fd(0)) < myEpsNull) {
1036     fd(0)=0.;
1037   }
1038
1039   k=n-1;
1040   for (i=1; i<k; i++) {
1041     fd(i)=.5*(f(i+1)-f(i-1))/(t(i)-t(i-1));
1042     if (fabs(fd(i)) < myEpsNull) {
1043       fd(i)=0.;
1044     }
1045   }
1046
1047   fd(n-1)=(f(n-1)-f(n-2))/(t(n-1)-t(n-2));
1048   if (fabs(fd(n-1)) < myEpsNull) {
1049     fd(n-1)=0.;
1050   }
1051   //
1052   // Finding the range where the derivatives have different signs
1053   // for neighbouring points
1054   for (i=1; i<n; i++) {
1055     Standard_Real fd1, fd2, t1, t2, fabsfd1, fabsfd2;
1056     Standard_Boolean bF1, bF2;
1057     t1 =t(i-1);
1058     t2 =t(i);
1059     fd1=fd(i-1);
1060     fd2=fd(i);
1061
1062     fabsfd1=fabs(fd1);
1063     bF1=fabsfd1 < myEpsNull;
1064     
1065     fabsfd2=fabs(fd2);
1066     bF2=fabsfd2 < myEpsNull;
1067     //
1068     if (fd1*fd2 < 0.) {
1069       tr=FindSimpleRoot(2, t1, t2, fd1);
1070       fr=DistanceFunction(tr);
1071       myPar1=tr;
1072       myParallel=Standard_False;
1073       break;
1074     }
1075     
1076     if (!bF1 && bF2) {
1077       tr=t2;
1078       fr=fd2;
1079       myPar1=tr;
1080       myParallel=Standard_False;
1081       break;
1082     }
1083     
1084     if (bF1 && !bF2) {
1085       tr=t1;
1086       fr=fd1;
1087       myPar1=tr;
1088       myParallel=Standard_False;
1089       break;
1090     }
1091
1092   }
1093 }
1094 //=======================================================================
1095 //function : RemoveIdenticalRoots 
1096 //purpose  : 
1097 //=======================================================================
1098   void IntTools_EdgeFace::RemoveIdenticalRoots()
1099 {
1100   Standard_Integer aNbRoots, j, k;
1101
1102   aNbRoots=mySequenceOfRoots.Length();
1103   for (j=1; j<=aNbRoots; j++) { 
1104     const IntTools_Root& aRj=mySequenceOfRoots(j);
1105     for (k=j+1; k<=aNbRoots; k++) {
1106       const IntTools_Root& aRk=mySequenceOfRoots(k);
1107       
1108       Standard_Real aTj, aTk, aDistance;
1109       gp_Pnt aPj, aPk;
1110
1111       aTj=aRj.Root();
1112       aTk=aRk.Root();
1113
1114       myC.D0(aTj, aPj);
1115       myC.D0(aTk, aPk);
1116
1117       aDistance=aPj.Distance(aPk);
1118       if (aDistance < myCriteria) {
1119         mySequenceOfRoots.Remove(k);
1120         aNbRoots=mySequenceOfRoots.Length();
1121       }
1122     }
1123   }
1124 }
1125
1126 //=======================================================================
1127 //function : CheckTouch 
1128 //purpose  : 
1129 //=======================================================================
1130   Standard_Boolean IntTools_EdgeFace::CheckTouch(const IntTools_CommonPrt& aCP,
1131                                                  Standard_Real&            aTx) 
1132 {
1133   Standard_Real aTF, aTL, Tol, U1f, U1l, V1f, V1l, af, al,aDist2, aMinDist2;
1134   Standard_Boolean theflag=Standard_False;
1135   Standard_Integer aNbExt, i, iLower ;
1136
1137   aCP.Range1(aTF, aTL);
1138
1139   //
1140   Standard_Real aCR;
1141   aCR=myC.Resolution(myCriteria);
1142   if((Abs(aTF - myRange.First()) < aCR) &&
1143      (Abs(aTL - myRange.Last())  < aCR)) {
1144     return theflag; // EDGE 
1145   }
1146   //
1147
1148   Tol = Precision::PConfusion();
1149
1150   const Handle(Geom_Curve)&  Curve   =BRep_Tool::Curve  (myC.Edge(), af, al);
1151   const Handle(Geom_Surface)& Surface=BRep_Tool::Surface(myS.Face());
1152   //   Surface->Bounds(U1f,U1l,V1f,V1l);
1153   U1f = myS.FirstUParameter();
1154   U1l = myS.LastUParameter();
1155   V1f = myS.FirstVParameter();
1156   V1l = myS.LastVParameter();
1157   
1158   GeomAdaptor_Curve   TheCurve   (Curve,aTF, aTL);
1159   GeomAdaptor_Surface TheSurface (Surface, U1f, U1l, V1f, V1l); 
1160                                  
1161   Extrema_ExtCS anExtrema (TheCurve, TheSurface, Tol, Tol);
1162
1163   aDist2 = 1.e100;
1164
1165   if(anExtrema.IsDone()) {
1166     aMinDist2 = aDist2;
1167
1168     if(!anExtrema.IsParallel()) {
1169       aNbExt=anExtrema.NbExt();
1170       
1171       if(aNbExt > 0) {
1172         iLower=1;
1173         for (i=1; i<=aNbExt; i++) {
1174           aDist2=anExtrema.SquareDistance(i);
1175           if (aDist2 < aMinDist2) {
1176             aMinDist2=aDist2;
1177             iLower=i;
1178           }
1179         }
1180         aDist2=anExtrema.SquareDistance(iLower);
1181         Extrema_POnCurv aPOnC;
1182         Extrema_POnSurf aPOnS;
1183         anExtrema.Points(iLower, aPOnC, aPOnS);
1184         aTx=aPOnC.Parameter();
1185       }
1186       else {
1187         // modified by NIZHNY-MKK  Thu Jul 21 11:35:32 2005.BEGIN
1188         IntCurveSurface_HInter anExactIntersector;
1189   
1190         Handle(GeomAdaptor_HCurve) aCurve     = new GeomAdaptor_HCurve(TheCurve);
1191         Handle(GeomAdaptor_HSurface) aSurface = new GeomAdaptor_HSurface(TheSurface);
1192         
1193         anExactIntersector.Perform(aCurve, aSurface);
1194
1195         if(anExactIntersector.IsDone()) {
1196           Standard_Integer i = 0;
1197
1198           for(i = 1; i <= anExactIntersector.NbPoints(); i++) {
1199             const IntCurveSurface_IntersectionPoint& aPoint = anExactIntersector.Point(i);
1200       
1201             if((aPoint.W() >= aTF) && (aPoint.W() <= aTL)) {
1202               aDist2=0.;
1203               aTx = aPoint.W();
1204             }
1205           }
1206         }
1207         // modified by NIZHNY-MKK  Thu Jul 21 11:35:40 2005.END
1208       }
1209     }
1210     else {
1211       return theflag;
1212     }
1213   }
1214
1215   Standard_Real aBoundaryDist;
1216
1217   aBoundaryDist = DistanceFunction(aTF) + myCriteria;
1218   if(aBoundaryDist * aBoundaryDist < aDist2) {
1219     aDist2 = aBoundaryDist * aBoundaryDist;
1220     aTx = aTF;
1221   }
1222   
1223   aBoundaryDist = DistanceFunction(aTL) + myCriteria;
1224   if(aBoundaryDist * aBoundaryDist < aDist2) {
1225     aDist2 = aBoundaryDist * aBoundaryDist;
1226     aTx = aTL;
1227   }
1228
1229   Standard_Real aParameter = (aTF + aTL) * 0.5;
1230   aBoundaryDist = DistanceFunction(aParameter) + myCriteria;
1231   if(aBoundaryDist * aBoundaryDist < aDist2) {
1232     aDist2 = aBoundaryDist * aBoundaryDist;
1233     aTx = aParameter;
1234   }
1235
1236   if(aDist2 > myCriteria * myCriteria) {
1237     return theflag;
1238   }
1239   
1240   if (fabs (aTx-aTF) < myEpsT) {
1241     return !theflag;
1242   }
1243
1244   if (fabs (aTx-aTL) < myEpsT) {
1245     return !theflag;
1246   }
1247
1248   if (aTx>aTF && aTx<aTL) {
1249     return !theflag;
1250   }
1251
1252   return theflag;
1253 }
1254
1255
1256 //=======================================================================
1257 //function : Perform
1258 //purpose  : 
1259 //=======================================================================
1260   void IntTools_EdgeFace::Perform() 
1261 {
1262   Standard_Integer i, aNb;
1263   IntTools_CommonPrt aCommonPrt;
1264   //
1265   aCommonPrt.SetEdge1(myEdge);
1266   //
1267   myErrorStatus=0;
1268   CheckData();
1269   if (myErrorStatus) {
1270     return;
1271   }
1272   //
1273   if (myContext.IsNull()) {
1274     myContext=new IntTools_Context;
1275   }
1276   //
1277   myIsDone = Standard_False;
1278   myC.Initialize(myEdge);
1279   GeomAbs_CurveType aCurveType;
1280   aCurveType=myC.GetType();
1281   //
1282   // Prepare myCriteria
1283   if (aCurveType==GeomAbs_BSplineCurve||
1284         aCurveType==GeomAbs_BezierCurve) {
1285     //--- 5112
1286     Standard_Real diff1 = (myTolE/myTolF);
1287     Standard_Real diff2 = (myTolF/myTolE);
1288     if( diff1 > 100 || diff2 > 100 ) {
1289       myCriteria = Max(myTolE,myTolF);
1290     }
1291     else //--- 5112
1292       myCriteria=1.5*myTolE+myTolF;
1293   }
1294   else {
1295     myCriteria=myTolE+myTolF;
1296   }
1297
1298   myTmin=myRange.First();
1299   myTmax=myRange.Last();
1300
1301   myS.Initialize (myFace,Standard_True);
1302
1303   if(myContext.IsNull()) {
1304     myFClass2d.Init(myFace, 1.e-6);
1305   }
1306
1307   IntTools_BeanFaceIntersector anIntersector(myC, myS, myTolE, myTolF);
1308   anIntersector.SetBeanParameters(myRange.First(), myRange.Last());
1309   //
1310   anIntersector.SetContext(myContext);
1311   //
1312   anIntersector.Perform();
1313
1314   if(!anIntersector.IsDone()) {
1315     return;
1316   }
1317
1318   for(Standard_Integer r = 1; r <= anIntersector.Result().Length(); r++) {
1319     const IntTools_Range& aRange = anIntersector.Result().Value(r);
1320     
1321     if(IsProjectable(IntTools_Tools::IntermediatePoint(aRange.First(), aRange.Last()))) {
1322       aCommonPrt.SetRange1(aRange.First(), aRange.Last());
1323       mySeqOfCommonPrts.Append(aCommonPrt);
1324     }
1325   }
1326
1327   aNb = mySeqOfCommonPrts.Length();
1328
1329   for (i=1; i<=aNb; i++) {
1330     IntTools_CommonPrt& aCP=mySeqOfCommonPrts.ChangeValue(i);
1331     //
1332     Standard_Real aTx1, aTx2;
1333     gp_Pnt aPx1, aPx2;
1334     //
1335     aCP.Range1(aTx1, aTx2);
1336     myC.D0(aTx1, aPx1);
1337     myC.D0(aTx2, aPx2);
1338     aCP.SetBoundingPoints(aPx1, aPx2);
1339     //
1340     MakeType (aCP); 
1341   }
1342   {
1343     // Line\Cylinder's Common Parts treatement
1344     GeomAbs_CurveType   aCType;
1345     GeomAbs_SurfaceType aSType;
1346     TopAbs_ShapeEnum aType;
1347     Standard_Boolean bIsTouch;
1348     Standard_Real aTx;
1349
1350     aCType=myC.GetType();
1351     aSType=myS.GetType();
1352     
1353     if (aCType==GeomAbs_Line && aSType==GeomAbs_Cylinder) {
1354       for (i=1; i<=aNb; i++) {
1355         IntTools_CommonPrt& aCP=mySeqOfCommonPrts(i);
1356         aType=aCP.Type();
1357         if (aType==TopAbs_EDGE) {
1358           bIsTouch=CheckTouch (aCP, aTx);
1359           if (bIsTouch) {
1360             aCP.SetType(TopAbs_VERTEX);
1361             aCP.SetVertexParameter1(aTx);
1362             aCP.SetRange1 (aTx, aTx);
1363           }
1364         }
1365         if (aType==TopAbs_VERTEX) {
1366           bIsTouch=CheckTouchVertex (aCP, aTx);
1367           if (bIsTouch) {
1368             aCP.SetVertexParameter1(aTx);
1369             aCP.SetRange1 (aTx, aTx);
1370           }
1371         }
1372       }
1373     }
1374
1375     // Circle\Plane's Common Parts treatement
1376     
1377     if (aCType==GeomAbs_Circle && aSType==GeomAbs_Plane) {
1378       Standard_Boolean bIsCoplanar, bIsRadius;
1379       bIsCoplanar=IsCoplanar(myC, myS);
1380       bIsRadius=IsRadius(myC, myS);
1381       if (!bIsCoplanar && !bIsRadius) {
1382         for (i=1; i<=aNb; i++) {
1383           IntTools_CommonPrt& aCP=mySeqOfCommonPrts(i);
1384           aType=aCP.Type();
1385           if (aType==TopAbs_EDGE) {
1386             bIsTouch=CheckTouch (aCP, aTx);
1387             if (bIsTouch) {
1388               aCP.SetType(TopAbs_VERTEX);
1389               aCP.SetVertexParameter1(aTx);
1390               aCP.SetRange1 (aTx, aTx);
1391             }
1392           }
1393         }
1394       }
1395     }
1396   }
1397   myIsDone=Standard_True;
1398 }
1399
1400 //
1401 // myErrorStatus
1402 // 1 - the method Perform() is not invoked  
1403 // 2,3,4,5 -the method CheckData() fails
1404 // 6 - PrepareArgs() problems
1405 // 7 - No Projectable ranges
1406 // 8,9 - PrepareArgs() problems occured inside  projectable Ranges
1407 // 11 - can't fill array  aFunc(i) in PrepareArgsFuncArrays 
1408
1409
1410 //=======================================================================
1411 //function : CheckTouch 
1412 //purpose  : 
1413 //=======================================================================
1414   Standard_Boolean IntTools_EdgeFace::CheckTouchVertex (const IntTools_CommonPrt& aCP,
1415                                                         Standard_Real& aTx) 
1416 {
1417   Standard_Real aTF, aTL, Tol, U1f,U1l,V1f,V1l, af, al,aDist2, aMinDist2, aTm, aDist2New;
1418   Standard_Real aEpsT;
1419   Standard_Boolean theflag=Standard_False;
1420   Standard_Integer aNbExt, i, iLower ;
1421
1422   aCP.Range1(aTF, aTL);
1423   aEpsT=8.e-5;
1424   aTm=0.5*(aTF+aTL);
1425   aDist2=DistanceFunction(aTm);
1426   aDist2 *= aDist2;
1427
1428   Tol = Precision::PConfusion();
1429
1430   const Handle(Geom_Curve)&  Curve   =BRep_Tool::Curve  (myC.Edge(), af, al);
1431   const Handle(Geom_Surface)& Surface=BRep_Tool::Surface(myS.Face());
1432
1433   Surface->Bounds(U1f,U1l,V1f,V1l);
1434   
1435   GeomAdaptor_Curve   TheCurve   (Curve,aTF, aTL);
1436   GeomAdaptor_Surface TheSurface (Surface, U1f, U1l, V1f, V1l); 
1437                                  
1438   Extrema_ExtCS anExtrema (TheCurve, TheSurface, Tol, Tol);
1439    
1440   if(!anExtrema.IsDone()) {
1441     return theflag;
1442   }
1443   if (anExtrema.IsParallel()) {
1444     return theflag;
1445   }
1446   
1447   aNbExt=anExtrema.NbExt() ;
1448   if (!aNbExt) {
1449      return theflag;
1450   }
1451
1452   iLower=1;
1453   aMinDist2=1.e100;
1454   for (i=1; i<=aNbExt; ++i) {
1455     aDist2=anExtrema.SquareDistance(i);
1456     if (aDist2 < aMinDist2) {
1457       aMinDist2=aDist2;
1458       iLower=i;
1459     }
1460   }
1461
1462   aDist2New=anExtrema.SquareDistance(iLower);
1463   
1464   if (aDist2New > aDist2) {
1465     aTx=aTm;
1466     return !theflag;
1467   }
1468   
1469   if (aDist2New > myCriteria * myCriteria) {
1470     return theflag;
1471   }
1472
1473   Extrema_POnCurv aPOnC;
1474   Extrema_POnSurf aPOnS;
1475   anExtrema.Points(iLower, aPOnC, aPOnS);
1476
1477  
1478   aTx=aPOnC.Parameter();
1479   ///
1480   if (fabs (aTx-aTF) < aEpsT) {
1481     return theflag;
1482   }
1483
1484   if (fabs (aTx-aTL) < aEpsT) {
1485     return theflag;
1486   }
1487
1488   if (aTx>aTF && aTx<aTL) {
1489     return !theflag;
1490   }
1491
1492   return theflag;
1493 }
1494
1495
1496 //=======================================================================
1497 //function :  IsCoplanar
1498 //purpose  : 
1499 //=======================================================================
1500 Standard_Boolean IsCoplanar (const BRepAdaptor_Curve& aCurve ,
1501                                const BRepAdaptor_Surface& aSurface)
1502 {
1503   Standard_Boolean bFlag=Standard_False;
1504
1505   GeomAbs_CurveType   aCType;
1506   GeomAbs_SurfaceType aSType;
1507
1508   aCType=aCurve.GetType();
1509   aSType=aSurface.GetType();
1510     
1511   if (aCType==GeomAbs_Circle && aSType==GeomAbs_Plane) {
1512     gp_Circ aCirc=aCurve.Circle();
1513     const gp_Ax1& anAx1=aCirc.Axis();
1514     const gp_Dir& aDirAx1=anAx1.Direction();
1515     
1516     gp_Pln aPln=aSurface.Plane();
1517     const gp_Ax1& anAx=aPln.Axis();
1518     const gp_Dir& aDirPln=anAx.Direction();
1519
1520     bFlag=IntTools_Tools::IsDirsCoinside(aDirAx1, aDirPln);
1521   }
1522   return bFlag;
1523 }
1524 //=======================================================================
1525 //function :  IsRadius
1526 //purpose  : 
1527 //=======================================================================
1528 Standard_Boolean IsRadius (const BRepAdaptor_Curve& aCurve ,
1529                            const BRepAdaptor_Surface& aSurface)
1530 {
1531   Standard_Boolean bFlag=Standard_False;
1532
1533   GeomAbs_CurveType   aCType;
1534   GeomAbs_SurfaceType aSType;
1535
1536   aCType=aCurve.GetType();
1537   aSType=aSurface.GetType();
1538     
1539   if (aCType==GeomAbs_Circle && aSType==GeomAbs_Plane) {
1540     gp_Circ aCirc=aCurve.Circle();
1541     const gp_Pnt aCenter=aCirc.Location();
1542     Standard_Real aR=aCirc.Radius();
1543     gp_Pln aPln=aSurface.Plane();
1544     Standard_Real aD=aPln.Distance(aCenter);
1545     if (fabs (aD-aR) < 1.e-7) {
1546       return !bFlag;
1547     }
1548   }
1549   return bFlag;
1550 }
1551 //
1552 //=======================================================================
1553 //function :  AdaptiveDiscret
1554 //purpose  : 
1555 //=======================================================================
1556 Standard_Integer AdaptiveDiscret (const Standard_Integer iDiscret,
1557                                   const BRepAdaptor_Curve& aCurve ,
1558                                   const BRepAdaptor_Surface& aSurface)
1559 {
1560   Standard_Integer iDiscretNew;
1561
1562   iDiscretNew=iDiscret;
1563
1564   GeomAbs_CurveType   aCType;
1565   GeomAbs_SurfaceType aSType;
1566
1567   aCType=aCurve.GetType();
1568   aSType=aSurface.GetType();
1569     
1570   if (aSType==GeomAbs_Cylinder) {
1571    Standard_Real aELength, aRadius, dL, dLR;
1572
1573    aELength=IntTools::Length(aCurve.Edge());
1574    dL=aELength/iDiscret;
1575    
1576    gp_Cylinder aCylinder=aSurface.Cylinder();
1577    aRadius=aCylinder.Radius();
1578    dLR=2*aRadius;
1579
1580    iDiscretNew=(Standard_Integer)(aELength/dLR);
1581    
1582    if (iDiscretNew<iDiscret) {
1583      iDiscretNew=iDiscret;
1584    }
1585      
1586   }
1587   return iDiscretNew;
1588 }
1589
1590
1591 #ifdef WNT
1592 #pragma warning ( default : 4101 )
1593 #endif