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