7952ffe24c1a9cbcc651d3a08fc49dbe995478bf
[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   return 0;
944 }
945
946
947 //=======================================================================
948 //function : IsIntersection
949 //purpose  : 
950 //=======================================================================
951   void IntTools_EdgeFace::IsIntersection (const Standard_Real ta, 
952                                           const Standard_Real tb) 
953 {
954   IntTools_CArray1OfReal anArgs, aFunc;
955   Standard_Integer i, aNb, pri, aCnt=0;
956   //
957   Standard_Integer aCntIncreasing=1, aCntDecreasing=1;
958   Standard_Real t, f, f1;
959   //
960   // Prepare values of arguments for the interval [ta, tb]
961   pri=IntTools::PrepareArgs (myC, tb, ta, myDiscret, myDeflection, anArgs);
962   aNb=anArgs.Length();
963   
964   aFunc.Resize(aNb);
965   for (i=0; i<aNb; i++) {
966     t=anArgs(i);
967     
968     f1=DistanceFunction(t);
969     f=f1+myCriteria; 
970
971     if (fabs(f1) < myEpsNull) { 
972       aCnt++;
973       f=0.;
974     }
975     aFunc(i)=f;
976     //
977     if (i) {
978       if (aFunc(i)>aFunc(i-1)) {
979         aCntIncreasing++;
980       }
981       if (aFunc(i)<aFunc(i-1)) {
982         aCntDecreasing++;
983       }
984     }
985     //
986   }
987
988   if (aCnt==aNb) {
989     myParallel=Standard_True;
990     return;
991   }
992   
993   FindDerivativeRoot(anArgs, aFunc);
994
995   //
996   if (myParallel) {
997     if (!(myC.GetType()==GeomAbs_Line 
998           && 
999           myS.GetType()==GeomAbs_Cylinder)) {
1000       if (aCntDecreasing==aNb) {
1001         myPar1=anArgs(aNb-1);
1002         myParallel=Standard_False;
1003       }
1004       if (aCntIncreasing==aNb) {
1005         myPar1=anArgs(0);
1006         myParallel=Standard_False;
1007       }
1008     }
1009   }
1010   //
1011   return ;
1012 }
1013
1014 //=======================================================================
1015 //function : FindDerivativeRoot
1016 //purpose  : 
1017 //=======================================================================
1018   void IntTools_EdgeFace::FindDerivativeRoot(const IntTools_CArray1OfReal& t,
1019                                              const IntTools_CArray1OfReal& f)  
1020 {
1021   Standard_Integer i, n, k;
1022   Standard_Real fr, tr;
1023   IntTools_CArray1OfReal fd;
1024   TColStd_SequenceOfReal aTSeq, aFSeq;  
1025   
1026   myPar1=0.;
1027   myParallel=Standard_True;
1028   
1029   n=t.Length();
1030   fd.Resize(n+1);
1031   //
1032   // Table of derivatives
1033   fd(0)=(f(1)-f(0))/(t(1)-t(0));
1034   if (fabs(fd(0)) < myEpsNull) {
1035     fd(0)=0.;
1036   }
1037
1038   k=n-1;
1039   for (i=1; i<k; i++) {
1040     fd(i)=.5*(f(i+1)-f(i-1))/(t(i)-t(i-1));
1041     if (fabs(fd(i)) < myEpsNull) {
1042       fd(i)=0.;
1043     }
1044   }
1045
1046   fd(n-1)=(f(n-1)-f(n-2))/(t(n-1)-t(n-2));
1047   if (fabs(fd(n-1)) < myEpsNull) {
1048     fd(n-1)=0.;
1049   }
1050   //
1051   // Finding the range where the derivatives have different signs
1052   // for neighbouring points
1053   for (i=1; i<n; i++) {
1054     Standard_Real fd1, fd2, t1, t2, fabsfd1, fabsfd2;
1055     Standard_Boolean bF1, bF2;
1056     t1 =t(i-1);
1057     t2 =t(i);
1058     fd1=fd(i-1);
1059     fd2=fd(i);
1060
1061     fabsfd1=fabs(fd1);
1062     bF1=fabsfd1 < myEpsNull;
1063     
1064     fabsfd2=fabs(fd2);
1065     bF2=fabsfd2 < myEpsNull;
1066     //
1067     if (fd1*fd2 < 0.) {
1068       tr=FindSimpleRoot(2, t1, t2, fd1);
1069       fr=DistanceFunction(tr);
1070       myPar1=tr;
1071       myParallel=Standard_False;
1072       break;
1073     }
1074     
1075     if (!bF1 && bF2) {
1076       tr=t2;
1077       fr=fd2;
1078       myPar1=tr;
1079       myParallel=Standard_False;
1080       break;
1081     }
1082     
1083     if (bF1 && !bF2) {
1084       tr=t1;
1085       fr=fd1;
1086       myPar1=tr;
1087       myParallel=Standard_False;
1088       break;
1089     }
1090
1091   }
1092 }
1093 //=======================================================================
1094 //function : RemoveIdenticalRoots 
1095 //purpose  : 
1096 //=======================================================================
1097   void IntTools_EdgeFace::RemoveIdenticalRoots()
1098 {
1099   Standard_Integer aNbRoots, j, k;
1100
1101   aNbRoots=mySequenceOfRoots.Length();
1102   for (j=1; j<=aNbRoots; j++) { 
1103     const IntTools_Root& aRj=mySequenceOfRoots(j);
1104     for (k=j+1; k<=aNbRoots; k++) {
1105       const IntTools_Root& aRk=mySequenceOfRoots(k);
1106       
1107       Standard_Real aTj, aTk, aDistance;
1108       gp_Pnt aPj, aPk;
1109
1110       aTj=aRj.Root();
1111       aTk=aRk.Root();
1112
1113       myC.D0(aTj, aPj);
1114       myC.D0(aTk, aPk);
1115
1116       aDistance=aPj.Distance(aPk);
1117       if (aDistance < myCriteria) {
1118         mySequenceOfRoots.Remove(k);
1119         aNbRoots=mySequenceOfRoots.Length();
1120       }
1121     }
1122   }
1123 }
1124
1125 //=======================================================================
1126 //function : CheckTouch 
1127 //purpose  : 
1128 //=======================================================================
1129   Standard_Boolean IntTools_EdgeFace::CheckTouch(const IntTools_CommonPrt& aCP,
1130                                                  Standard_Real&            aTx) 
1131 {
1132   Standard_Real aTF, aTL, Tol, U1f, U1l, V1f, V1l, af, al,aDist2, aMinDist2;
1133   Standard_Boolean theflag=Standard_False;
1134   Standard_Integer aNbExt, i, iLower ;
1135
1136   aCP.Range1(aTF, aTL);
1137
1138   //
1139   Standard_Real aCR;
1140   aCR=myC.Resolution(myCriteria);
1141   if((Abs(aTF - myRange.First()) < aCR) &&
1142      (Abs(aTL - myRange.Last())  < aCR)) {
1143     return theflag; // EDGE 
1144   }
1145   //
1146
1147   Tol = Precision::PConfusion();
1148
1149   const Handle(Geom_Curve)&  Curve   =BRep_Tool::Curve  (myC.Edge(), af, al);
1150   const Handle(Geom_Surface)& Surface=BRep_Tool::Surface(myS.Face());
1151   //   Surface->Bounds(U1f,U1l,V1f,V1l);
1152   U1f = myS.FirstUParameter();
1153   U1l = myS.LastUParameter();
1154   V1f = myS.FirstVParameter();
1155   V1l = myS.LastVParameter();
1156   
1157   GeomAdaptor_Curve   TheCurve   (Curve,aTF, aTL);
1158   GeomAdaptor_Surface TheSurface (Surface, U1f, U1l, V1f, V1l); 
1159                                  
1160   Extrema_ExtCS anExtrema (TheCurve, TheSurface, Tol, Tol);
1161
1162   aDist2 = 1.e100;
1163
1164   if(anExtrema.IsDone()) {
1165     aMinDist2 = aDist2;
1166
1167     if(!anExtrema.IsParallel()) {
1168       aNbExt=anExtrema.NbExt();
1169       
1170       if(aNbExt > 0) {
1171         iLower=1;
1172         for (i=1; i<=aNbExt; i++) {
1173           aDist2=anExtrema.SquareDistance(i);
1174           if (aDist2 < aMinDist2) {
1175             aMinDist2=aDist2;
1176             iLower=i;
1177           }
1178         }
1179         aDist2=anExtrema.SquareDistance(iLower);
1180         Extrema_POnCurv aPOnC;
1181         Extrema_POnSurf aPOnS;
1182         anExtrema.Points(iLower, aPOnC, aPOnS);
1183         aTx=aPOnC.Parameter();
1184       }
1185       else {
1186         // modified by NIZHNY-MKK  Thu Jul 21 11:35:32 2005.BEGIN
1187         IntCurveSurface_HInter anExactIntersector;
1188   
1189         Handle(GeomAdaptor_HCurve) aCurve     = new GeomAdaptor_HCurve(TheCurve);
1190         Handle(GeomAdaptor_HSurface) aSurface = new GeomAdaptor_HSurface(TheSurface);
1191         
1192         anExactIntersector.Perform(aCurve, aSurface);
1193
1194         if(anExactIntersector.IsDone()) {
1195           Standard_Integer i = 0;
1196
1197           for(i = 1; i <= anExactIntersector.NbPoints(); i++) {
1198             const IntCurveSurface_IntersectionPoint& aPoint = anExactIntersector.Point(i);
1199       
1200             if((aPoint.W() >= aTF) && (aPoint.W() <= aTL)) {
1201               aDist2=0.;
1202               aTx = aPoint.W();
1203             }
1204           }
1205         }
1206         // modified by NIZHNY-MKK  Thu Jul 21 11:35:40 2005.END
1207       }
1208     }
1209     else {
1210       return theflag;
1211     }
1212   }
1213
1214   Standard_Real aBoundaryDist;
1215
1216   aBoundaryDist = DistanceFunction(aTF) + myCriteria;
1217   if(aBoundaryDist * aBoundaryDist < aDist2) {
1218     aDist2 = aBoundaryDist * aBoundaryDist;
1219     aTx = aTF;
1220   }
1221   
1222   aBoundaryDist = DistanceFunction(aTL) + myCriteria;
1223   if(aBoundaryDist * aBoundaryDist < aDist2) {
1224     aDist2 = aBoundaryDist * aBoundaryDist;
1225     aTx = aTL;
1226   }
1227
1228   Standard_Real aParameter = (aTF + aTL) * 0.5;
1229   aBoundaryDist = DistanceFunction(aParameter) + myCriteria;
1230   if(aBoundaryDist * aBoundaryDist < aDist2) {
1231     aDist2 = aBoundaryDist * aBoundaryDist;
1232     aTx = aParameter;
1233   }
1234
1235   if(aDist2 > myCriteria * myCriteria) {
1236     return theflag;
1237   }
1238   
1239   if (fabs (aTx-aTF) < myEpsT) {
1240     return !theflag;
1241   }
1242
1243   if (fabs (aTx-aTL) < myEpsT) {
1244     return !theflag;
1245   }
1246
1247   if (aTx>aTF && aTx<aTL) {
1248     return !theflag;
1249   }
1250
1251   return theflag;
1252 }
1253
1254
1255 //=======================================================================
1256 //function : Perform
1257 //purpose  : 
1258 //=======================================================================
1259   void IntTools_EdgeFace::Perform() 
1260 {
1261   Standard_Integer i, aNb;
1262   IntTools_CommonPrt aCommonPrt;
1263   //
1264   aCommonPrt.SetEdge1(myEdge);
1265   //
1266   myErrorStatus=0;
1267   CheckData();
1268   if (myErrorStatus) {
1269     return;
1270   }
1271   //
1272   if (myContext.IsNull()) {
1273     myContext=new IntTools_Context;
1274   }
1275   //
1276   myIsDone = Standard_False;
1277   myC.Initialize(myEdge);
1278   GeomAbs_CurveType aCurveType;
1279   aCurveType=myC.GetType();
1280   //
1281   // Prepare myCriteria
1282   if (aCurveType==GeomAbs_BSplineCurve||
1283         aCurveType==GeomAbs_BezierCurve) {
1284     //--- 5112
1285     Standard_Real diff1 = (myTolE/myTolF);
1286     Standard_Real diff2 = (myTolF/myTolE);
1287     if( diff1 > 100 || diff2 > 100 ) {
1288       myCriteria = Max(myTolE,myTolF);
1289     }
1290     else //--- 5112
1291       myCriteria=1.5*myTolE+myTolF;
1292   }
1293   else {
1294     myCriteria=myTolE+myTolF;
1295   }
1296
1297   myTmin=myRange.First();
1298   myTmax=myRange.Last();
1299
1300   myS.Initialize (myFace,Standard_True);
1301
1302   if(myContext.IsNull()) {
1303     myFClass2d.Init(myFace, 1.e-6);
1304   }
1305
1306   IntTools_BeanFaceIntersector anIntersector(myC, myS, myTolE, myTolF);
1307   anIntersector.SetBeanParameters(myRange.First(), myRange.Last());
1308   //
1309   anIntersector.SetContext(myContext);
1310   //
1311   anIntersector.Perform();
1312
1313   if(!anIntersector.IsDone()) {
1314     return;
1315   }
1316
1317   for(Standard_Integer r = 1; r <= anIntersector.Result().Length(); r++) {
1318     const IntTools_Range& aRange = anIntersector.Result().Value(r);
1319     
1320     if(IsProjectable(IntTools_Tools::IntermediatePoint(aRange.First(), aRange.Last()))) {
1321       aCommonPrt.SetRange1(aRange.First(), aRange.Last());
1322       mySeqOfCommonPrts.Append(aCommonPrt);
1323     }
1324   }
1325
1326   aNb = mySeqOfCommonPrts.Length();
1327
1328   for (i=1; i<=aNb; i++) {
1329     IntTools_CommonPrt& aCP=mySeqOfCommonPrts.ChangeValue(i);
1330     //
1331     Standard_Real aTx1, aTx2;
1332     gp_Pnt aPx1, aPx2;
1333     //
1334     aCP.Range1(aTx1, aTx2);
1335     myC.D0(aTx1, aPx1);
1336     myC.D0(aTx2, aPx2);
1337     aCP.SetBoundingPoints(aPx1, aPx2);
1338     //
1339     MakeType (aCP); 
1340   }
1341   {
1342     // Line\Cylinder's Common Parts treatement
1343     GeomAbs_CurveType   aCType;
1344     GeomAbs_SurfaceType aSType;
1345     TopAbs_ShapeEnum aType;
1346     Standard_Boolean bIsTouch;
1347     Standard_Real aTx;
1348
1349     aCType=myC.GetType();
1350     aSType=myS.GetType();
1351     
1352     if (aCType==GeomAbs_Line && aSType==GeomAbs_Cylinder) {
1353       for (i=1; i<=aNb; i++) {
1354         IntTools_CommonPrt& aCP=mySeqOfCommonPrts(i);
1355         aType=aCP.Type();
1356         if (aType==TopAbs_EDGE) {
1357           bIsTouch=CheckTouch (aCP, aTx);
1358           if (bIsTouch) {
1359             aCP.SetType(TopAbs_VERTEX);
1360             aCP.SetVertexParameter1(aTx);
1361             aCP.SetRange1 (aTx, aTx);
1362           }
1363         }
1364         if (aType==TopAbs_VERTEX) {
1365           bIsTouch=CheckTouchVertex (aCP, aTx);
1366           if (bIsTouch) {
1367             aCP.SetVertexParameter1(aTx);
1368             aCP.SetRange1 (aTx, aTx);
1369           }
1370         }
1371       }
1372     }
1373
1374     // Circle\Plane's Common Parts treatement
1375     
1376     if (aCType==GeomAbs_Circle && aSType==GeomAbs_Plane) {
1377       Standard_Boolean bIsCoplanar, bIsRadius;
1378       bIsCoplanar=IsCoplanar(myC, myS);
1379       bIsRadius=IsRadius(myC, myS);
1380       if (!bIsCoplanar && !bIsRadius) {
1381         for (i=1; i<=aNb; i++) {
1382           IntTools_CommonPrt& aCP=mySeqOfCommonPrts(i);
1383           aType=aCP.Type();
1384           if (aType==TopAbs_EDGE) {
1385             bIsTouch=CheckTouch (aCP, aTx);
1386             if (bIsTouch) {
1387               aCP.SetType(TopAbs_VERTEX);
1388               aCP.SetVertexParameter1(aTx);
1389               aCP.SetRange1 (aTx, aTx);
1390             }
1391           }
1392         }
1393       }
1394     }
1395   }
1396   myIsDone=Standard_True;
1397 }
1398
1399 //
1400 // myErrorStatus
1401 // 1 - the method Perform() is not invoked  
1402 // 2,3,4,5 -the method CheckData() fails
1403 // 6 - PrepareArgs() problems
1404 // 7 - No Projectable ranges
1405 // 8,9 - PrepareArgs() problems occured inside  projectable Ranges
1406 // 11 - can't fill array  aFunc(i) in PrepareArgsFuncArrays 
1407
1408
1409 //=======================================================================
1410 //function : CheckTouch 
1411 //purpose  : 
1412 //=======================================================================
1413   Standard_Boolean IntTools_EdgeFace::CheckTouchVertex (const IntTools_CommonPrt& aCP,
1414                                                         Standard_Real& aTx) 
1415 {
1416   Standard_Real aTF, aTL, Tol, U1f,U1l,V1f,V1l, af, al,aDist2, aMinDist2, aTm, aDist2New;
1417   Standard_Boolean theflag=Standard_False;
1418   Standard_Integer aNbExt, i, iLower ;
1419
1420   aCP.Range1(aTF, aTL);
1421
1422   aTm=0.5*(aTF+aTL);
1423   aDist2=DistanceFunction(aTm);
1424   aDist2 *= aDist2;
1425
1426   Tol = Precision::PConfusion();
1427
1428   const Handle(Geom_Curve)&  Curve   =BRep_Tool::Curve  (myC.Edge(), af, al);
1429   const Handle(Geom_Surface)& Surface=BRep_Tool::Surface(myS.Face());
1430
1431   Surface->Bounds(U1f,U1l,V1f,V1l);
1432   
1433   GeomAdaptor_Curve   TheCurve   (Curve,aTF, aTL);
1434   GeomAdaptor_Surface TheSurface (Surface, U1f, U1l, V1f, V1l); 
1435                                  
1436   Extrema_ExtCS anExtrema (TheCurve, TheSurface, Tol, Tol);
1437    
1438   if(!anExtrema.IsDone()) {
1439     return theflag;
1440   }
1441   if (anExtrema.IsParallel()) {
1442     return theflag;
1443   }
1444   
1445   aNbExt=anExtrema.NbExt() ;
1446   if (!aNbExt) {
1447      return theflag;
1448   }
1449
1450   iLower=1;
1451   aMinDist2=1.e100;
1452   for (i=1; i<=aNbExt; ++i) {
1453     aDist2=anExtrema.SquareDistance(i);
1454     if (aDist2 < aMinDist2) {
1455       aMinDist2=aDist2;
1456       iLower=i;
1457     }
1458   }
1459
1460   aDist2New=anExtrema.SquareDistance(iLower);
1461   
1462   if (aDist2New > aDist2) {
1463     aTx=aTm;
1464     return !theflag;
1465   }
1466   
1467   if (aDist2New > myCriteria * myCriteria) {
1468     return theflag;
1469   }
1470
1471   Extrema_POnCurv aPOnC;
1472   Extrema_POnSurf aPOnS;
1473   anExtrema.Points(iLower, aPOnC, aPOnS);
1474
1475   
1476   aTx=aPOnC.Parameter();
1477   
1478   if (fabs (aTx-aTF) < myEpsT) {
1479     return !theflag;
1480   }
1481
1482   if (fabs (aTx-aTL) < myEpsT) {
1483     return !theflag;
1484   }
1485
1486   if (aTx>aTF && aTx<aTL) {
1487     return !theflag;
1488   }
1489
1490   return theflag;
1491 }
1492
1493
1494 //=======================================================================
1495 //function :  IsCoplanar
1496 //purpose  : 
1497 //=======================================================================
1498 Standard_Boolean IsCoplanar (const BRepAdaptor_Curve& aCurve ,
1499                                const BRepAdaptor_Surface& aSurface)
1500 {
1501   Standard_Boolean bFlag=Standard_False;
1502
1503   GeomAbs_CurveType   aCType;
1504   GeomAbs_SurfaceType aSType;
1505
1506   aCType=aCurve.GetType();
1507   aSType=aSurface.GetType();
1508     
1509   if (aCType==GeomAbs_Circle && aSType==GeomAbs_Plane) {
1510     gp_Circ aCirc=aCurve.Circle();
1511     const gp_Ax1& anAx1=aCirc.Axis();
1512     const gp_Dir& aDirAx1=anAx1.Direction();
1513     
1514     gp_Pln aPln=aSurface.Plane();
1515     const gp_Ax1& anAx=aPln.Axis();
1516     const gp_Dir& aDirPln=anAx.Direction();
1517
1518     bFlag=IntTools_Tools::IsDirsCoinside(aDirAx1, aDirPln);
1519   }
1520   return bFlag;
1521 }
1522 //=======================================================================
1523 //function :  IsRadius
1524 //purpose  : 
1525 //=======================================================================
1526 Standard_Boolean IsRadius (const BRepAdaptor_Curve& aCurve ,
1527                            const BRepAdaptor_Surface& aSurface)
1528 {
1529   Standard_Boolean bFlag=Standard_False;
1530
1531   GeomAbs_CurveType   aCType;
1532   GeomAbs_SurfaceType aSType;
1533
1534   aCType=aCurve.GetType();
1535   aSType=aSurface.GetType();
1536     
1537   if (aCType==GeomAbs_Circle && aSType==GeomAbs_Plane) {
1538     gp_Circ aCirc=aCurve.Circle();
1539     const gp_Pnt aCenter=aCirc.Location();
1540     Standard_Real aR=aCirc.Radius();
1541     gp_Pln aPln=aSurface.Plane();
1542     Standard_Real aD=aPln.Distance(aCenter);
1543     if (fabs (aD-aR) < 1.e-7) {
1544       return !bFlag;
1545     }
1546   }
1547   return bFlag;
1548 }
1549 //
1550 //=======================================================================
1551 //function :  AdaptiveDiscret
1552 //purpose  : 
1553 //=======================================================================
1554 Standard_Integer AdaptiveDiscret (const Standard_Integer iDiscret,
1555                                   const BRepAdaptor_Curve& aCurve ,
1556                                   const BRepAdaptor_Surface& aSurface)
1557 {
1558   Standard_Integer iDiscretNew;
1559
1560   iDiscretNew=iDiscret;
1561
1562   GeomAbs_CurveType   aCType;
1563   GeomAbs_SurfaceType aSType;
1564
1565   aCType=aCurve.GetType();
1566   aSType=aSurface.GetType();
1567     
1568   if (aSType==GeomAbs_Cylinder) {
1569    Standard_Real aELength, aRadius, dL, dLR;
1570
1571    aELength=IntTools::Length(aCurve.Edge());
1572    dL=aELength/iDiscret;
1573    
1574    gp_Cylinder aCylinder=aSurface.Cylinder();
1575    aRadius=aCylinder.Radius();
1576    dLR=2*aRadius;
1577
1578    iDiscretNew=(Standard_Integer)(aELength/dLR);
1579    
1580    if (iDiscretNew<iDiscret) {
1581      iDiscretNew=iDiscret;
1582    }
1583      
1584   }
1585   return iDiscretNew;
1586 }
1587
1588
1589 #ifdef WNT
1590 #pragma warning ( default : 4101 )
1591 #endif