0025246: Make methods Intervals and NbIntervals const in Adaptor3d_Curve and it desce...
[occt.git] / src / Adaptor3d / Adaptor3d_CurveOnSurface.cxx
1 // Copyright (c) 1999-2014 OPEN CASCADE SAS
2 //
3 // This file is part of Open CASCADE Technology software library.
4 //
5 // This library is free software; you can redistribute it and/or modify it under
6 // the terms of the GNU Lesser General Public License version 2.1 as published
7 // by the Free Software Foundation, with special exception defined in the file
8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
9 // distribution for complete text of the license and disclaimer of any warranty.
10 //
11 // Alternatively, this file may be used under the terms of Open CASCADE
12 // commercial license or contractual agreement.
13
14 #define No_Standard_OutOfRange
15
16 #include <Adaptor3d_CurveOnSurface.ixx>
17
18 #include <Adaptor3d_HCurveOnSurface.hxx>
19 #include <gp_Pnt2d.hxx>
20 #include <gp_Vec2d.hxx>
21 #include <gp_Ax22d.hxx>
22 #include <gp_Lin2d.hxx>
23 #include <gp_Circ2d.hxx>
24 #include <gp_Elips2d.hxx>
25 #include <gp_Hypr2d.hxx>
26 #include <gp_Parab2d.hxx>
27 #include <Geom_BSplineSurface.hxx>
28 #include <Geom_SurfaceOfRevolution.hxx>
29 #include <Geom_SurfaceOfLinearExtrusion.hxx>
30 #include <Geom_OffsetSurface.hxx>
31 #include <Geom2d_BezierCurve.hxx>
32 #include <Geom2d_BSplineCurve.hxx>
33 #include <Precision.hxx>
34 #include <Standard_Assert.hxx>
35 #include <TColgp_Array1OfPnt2d.hxx>
36 #include <TColgp_Array1OfPnt.hxx>
37 #include <TColStd_Array1OfReal.hxx>
38 #include <TColStd_Array1OfInteger.hxx>
39 #include <TColStd_HSequenceOfReal.hxx>
40 #include <Standard_NotImplemented.hxx>
41 #include <ElCLib.hxx>
42 #include <ElSLib.hxx>
43 #include <Adaptor3d_InterFunc.hxx>
44 #include <math_FunctionRoots.hxx>
45 #include <ElSLib.hxx>
46
47 static gp_Pnt to3d(const gp_Pln& Pl, const gp_Pnt2d& P)
48 {
49   return ElSLib::Value(P.X(),P.Y(),Pl);
50 }
51
52 static gp_Vec to3d(const gp_Pln& Pl, const gp_Vec2d& V)
53 {
54   gp_Vec Vx = Pl.XAxis().Direction();
55   gp_Vec Vy = Pl.YAxis().Direction();
56   Vx.Multiply(V.X());
57   Vy.Multiply(V.Y());
58   Vx.Add(Vy);
59   return Vx;
60 }
61
62 static gp_Ax2 to3d(const gp_Pln& Pl, const gp_Ax22d& A)
63 {
64   gp_Pnt P = to3d(Pl,A.Location());
65   gp_Vec VX = to3d(Pl,A.XAxis().Direction());
66   gp_Vec VY = to3d(Pl,A.YAxis().Direction());
67   return gp_Ax2(P,VX.Crossed(VY),VX);
68 }
69
70 static gp_Circ to3d(const gp_Pln& Pl, const gp_Circ2d& C)
71 {
72   return gp_Circ(to3d(Pl,C.Axis()),C.Radius());
73 }
74
75 static gp_Elips to3d(const gp_Pln& Pl, const gp_Elips2d& E)
76 {
77   return gp_Elips(to3d(Pl,E.Axis()),E.MajorRadius(),E.MinorRadius());
78 }
79
80 static gp_Hypr to3d(const gp_Pln& Pl, const gp_Hypr2d& H)
81 {
82   return gp_Hypr(to3d(Pl,H.Axis()),H.MajorRadius(),H.MinorRadius());
83 }
84
85 static gp_Parab to3d(const gp_Pln& Pl, const gp_Parab2d& P)
86 {
87   return gp_Parab(to3d(Pl,P.Axis()),P.Focal());
88 }
89
90 static gp_Vec SetLinearForm(const gp_Vec2d DW, const gp_Vec2d D2W,const gp_Vec2d D3W,
91                             const gp_Vec D1U,  const gp_Vec D1V,  const gp_Vec D2U,
92                             const gp_Vec D2V,  const gp_Vec D2UV, const gp_Vec D3U,
93                             const gp_Vec D3V,  const gp_Vec D3UUV,const gp_Vec D3UVV)
94 {gp_Vec V31, V32, V33, V34,V3 ;
95  V31.SetLinearForm(DW.X(),D1U,
96                    D2W.X()*DW.X(),D2U,
97                    D2W.X()*DW.Y(),D2UV);
98  V31.SetLinearForm(D3W.Y(),D1V,
99                    D2W.Y()*DW.X(),D2UV,
100                    D2W.Y()*DW.Y(),D2V,
101                    V31);
102  V32.SetLinearForm(DW.X()*DW.X()*DW.Y(),D3UUV,
103                    DW.X()*DW.Y()*DW.Y(),D3UVV);
104  V32.SetLinearForm(D2W.X()*DW.Y()+DW.X()*D2W.Y(),D2UV,
105                    DW.X()*DW.Y()*DW.Y(),D3UVV,
106                    V32);
107  V33.SetLinearForm(2*D2W.X()*DW.X(),D2U,
108                    DW.X()*DW.X()*DW.X(),D3U,
109                    DW.X()*DW.X()*DW.Y(),D3UUV);
110  
111  V34.SetLinearForm(2*D2W.Y()*DW.Y(),D2V,
112                    DW.Y()*DW.Y()*DW.X(),D3UVV,
113                    DW.Y()*DW.Y()*DW.Y(),D3V);
114  V3.SetLinearForm(1,V31,2,V32,1,V33,V34);
115  return V3;
116 }
117
118 //=======================================================================
119 static void CompareBounds(gp_Pnt2d& P1,
120                           gp_Pnt2d& P2)//SVV
121 {
122    Standard_Real Lx = P1.X(),Ly = P1.Y();
123    Standard_Real Rx = P2.X(),Ry = P2.Y();
124    
125    if (Lx > Rx) { P1.SetX(Rx); P2.SetX(Lx);}
126    if (Ly > Ry) { P1.SetY(Ry); P2.SetY(Ly);}
127 }
128
129 //=======================================================================
130 //function :Hunt
131 //purpose  : 
132 //=======================================================================
133 static void Hunt(const TColStd_Array1OfReal& Arr,
134                  const Standard_Real Coord,
135                  Standard_Integer& Iloc)
136 {//Warning: Hunt is used to find number of knot which equals co-ordinate component,
137   //        when co-ordinate component definitly equals a knot only.
138   Standard_Real Tol=Precision::PConfusion()/10;
139   Standard_Integer i=1; 
140   while((i <= Arr.Upper()) && (Abs(Coord - Arr(i)) > Tol)){
141     i++;}
142  
143   if(Abs(Coord - Arr(i)) < Tol)
144     Iloc = i;
145   else
146     if(Abs(Coord - Arr(i)) > Tol) 
147       Standard_NotImplemented::Raise("Adaptor3d_CurveOnSurface:Hunt");
148 }
149
150 //=======================================================================
151 //function :ReverseParam
152 //purpose  : 
153 //=======================================================================
154
155 static void ReverseParam(const Standard_Real In1,
156                          const Standard_Real In2,
157                          Standard_Real& Out1,
158                          Standard_Real& Out2 )
159 {
160
161   if(In1>In2)  {Out1=In2;
162                 Out2=In1;}
163   else         {Out1=In1;
164                 Out2=In2;}
165 }
166 //=======================================================================
167 //function :ReverseParam
168 //purpose  : 
169 //=======================================================================
170
171 static void ReverseParam(const Standard_Integer In1,
172                          const Standard_Integer In2,
173                          Standard_Integer& Out1,
174                          Standard_Integer& Out2 )
175 {
176   if(In1>In2)  {Out1=In2;
177                 Out2=In1;}
178   else         {Out1=In1;
179                 Out2=In2;}
180 }
181
182 //=======================================================================
183 //function :FindBounds
184 //purpose  : 
185 //=======================================================================
186 static void FindBounds(const TColStd_Array1OfReal& Arr,
187                        const Standard_Real Coord,
188                        const Standard_Real Der,
189                        Standard_Integer& Bound1,
190                        Standard_Integer& Bound2,
191                        Standard_Boolean& DerNull)
192
193 {
194   Standard_Integer N=0;
195   Standard_Real Tol=Precision::PConfusion()/10;
196   Hunt(Arr,Coord,N);
197   DerNull=Standard_False;
198   
199   if(N==Bound1){ if(Abs(Der) > Tol) DerNull = Standard_False;
200                       if(Abs(Der)<= Tol) DerNull = Standard_True;
201                       Bound1=N;Bound2=N+1; return;
202                     }
203   if(N==Bound2){ if( Abs(Der) > Tol ) DerNull = Standard_False;
204                       if( Abs(Der)<= Tol ) DerNull = Standard_True;
205                       Bound1=N-1;Bound2=N; return;
206                     }
207   if((N!=Bound1)&&(N!=Bound2))  { 
208     if(Abs(Der) > Tol ) { 
209       if(Der>0) {Bound1=N;Bound2= N+1;}
210       else
211       if(Der<0){Bound1=N-1;Bound2=N;}
212       DerNull = Standard_False;
213     }
214     if(Abs(Der) <=Tol ) {
215       DerNull = Standard_True;
216       Bound1=N-1;
217       Bound2=N+1;
218     }
219   } 
220 }
221
222 //=======================================================================
223 //function :Locate1Coord
224 //purpose  : along BSpline curve
225 //=======================================================================
226
227 static void Locate1Coord(const Standard_Integer Index,
228                          const gp_Pnt2d& UV, 
229                          const gp_Vec2d& DUV,
230                          const Handle(Geom_BSplineCurve)& BSplC,
231                          gp_Pnt2d& LeftBot, gp_Pnt2d& RightTop)
232 {
233   Standard_Real Comp1=0, DComp1=0, cur, f = 0.0, l = 0.0;
234   Standard_Real Tol = Precision::PConfusion()/10;
235   Standard_Integer i = 1, Bnd1, Bnd2;
236   Standard_Boolean DIsNull= Standard_False; 
237   TColStd_Array1OfReal Arr(1,BSplC->NbKnots());  BSplC->Knots(Arr);
238   
239   if(Index==1) {    Comp1=UV.X();  DComp1=DUV.X(); } 
240   if(Index==2) {    Comp1=UV.Y();  DComp1=DUV.Y(); } 
241
242   Standard_Integer Lo = BSplC->FirstUKnotIndex(), Up = BSplC->LastUKnotIndex();
243
244   i = Lo;
245   while ( ( Abs(BSplC->Knot(i)-Comp1)>Tol )&& (i!=Up ) ) i++;
246   cur=BSplC->Knot(i);
247
248   if( Abs(Comp1-cur)<=Tol) {
249  
250     Bnd1 = Lo; Bnd2 = Up;
251     FindBounds(Arr,cur,DComp1,Bnd1,Bnd2,DIsNull); 
252     ReverseParam(Bnd1,Bnd2,Bnd1,Bnd2);
253     
254     if(DIsNull==Standard_False){
255       if(Index==1)  {LeftBot.SetX(BSplC->Knot(Bnd1));
256                      RightTop.SetX(BSplC->Knot(Bnd2));}
257       else 
258         if(Index==2){ LeftBot.SetY(BSplC->Knot(Bnd1));
259                       RightTop.SetY(BSplC->Knot(Bnd2)); }
260     }
261     else
262       if(DIsNull==Standard_True){
263         if( Abs( Comp1-(f=BSplC->Knot(Lo))) <= Tol)
264           {       
265             if(Index==1)   { LeftBot.SetX(BSplC->Knot(Lo));
266                              RightTop.SetX(BSplC->Knot(Lo+1));}
267             else if(Index==2)  { LeftBot.SetY(BSplC->Knot(Lo));
268                                  RightTop.SetY(BSplC->Knot(Lo+1));}
269            } else
270              if( Abs( Comp1-(l=BSplC->Knot(Up))) <= Tol)
271                {
272                  if(Index==1)   {  LeftBot.SetX(BSplC->Knot(Up-1));
273                                    RightTop.SetX(BSplC->Knot(Up));}
274                  else if(Index==2)  {LeftBot.SetY(BSplC->Knot(Up-1));
275                                      RightTop.SetY(BSplC->Knot(Up));}
276                }else
277                  if(Index==1)   {  LeftBot.SetX(BSplC->Knot(Bnd1));
278                                    RightTop.SetX(BSplC->Knot(Bnd2));}
279                  else if(Index==2)  {LeftBot.SetY(BSplC->Knot(Bnd1));
280                                      RightTop.SetY(BSplC->Knot(Bnd2));}
281       }
282   } 
283   else//*********if Coord != Knot
284     {
285       i=Lo;
286       while (i < Up) {
287         //if((f=BSplC->Knot(i))<Comp1 && (l=BSplC->Knot(i+1))>Comp1) break;
288         //skl 28.03.2002 for OCC233
289         f=BSplC->Knot(i);
290         l=BSplC->Knot(i+1);
291         if(f<Comp1 && l>Comp1) break;
292         i++; 
293       }
294       ReverseParam(f,l,f,l);
295      
296       if(i!=Up)   {
297         if(Abs(DComp1)<Tol) 
298           { if(Index==1) {LeftBot.SetX(f); RightTop.SetX(l);}else
299               if(Index==2) {LeftBot.SetY(f); RightTop.SetY(l); }
300           }else
301             if(Abs(DComp1)>Tol) 
302               { 
303                 if(Index==1) {
304                   if(DComp1>0) {LeftBot.SetX(Comp1); RightTop.SetX(l);}   else
305                     if(DComp1<0) {LeftBot.SetX(f); RightTop.SetX(Comp1);}
306                 }
307                 else
308                   if(Index==2) {
309                     if(DComp1>0) {LeftBot.SetY(Comp1); RightTop.SetY(l);} else
310                       if(DComp1<0) {LeftBot.SetY(f); RightTop.SetY(Comp1);};
311                   }
312               }
313       }else
314         if(i==Up)   {
315           if(Index==1) {LeftBot.SetX(Comp1); RightTop.SetX(BSplC->Knot(i));}else
316             if(Index==2) {LeftBot.SetY(Comp1); RightTop.SetY(BSplC->Knot(i)); }
317         }
318     }
319 }
320
321
322 //=======================================================================
323 //function :Locate1Coord
324 //purpose  : 
325 //=======================================================================
326
327 static void Locate1Coord(const Standard_Integer Index,
328         const gp_Pnt2d& UV, 
329         const gp_Vec2d& DUV,
330         const Handle(Geom_BSplineSurface)& BSplS,
331         Standard_Boolean& DIsNull,
332         gp_Pnt2d& LeftBot, 
333         gp_Pnt2d& RightTop)
334 {
335   Standard_Real Comp1=0,DComp1=0; 
336   Standard_Real Tol = Precision::PConfusion()/10;
337   Standard_Integer i=1, Up=0, Up1, Up2, Down=0, Down1, Down2;
338   Standard_Real cur = 0.;
339
340   DIsNull= Standard_False; 
341
342   Up1 =  BSplS->LastUKnotIndex();
343   Down1 = BSplS->FirstUKnotIndex();
344   Up2 =  BSplS->LastVKnotIndex();
345   Down2 = BSplS->FirstVKnotIndex();
346
347   if(Index==1)
348   { 
349     i = Down1; 
350     Comp1 =  UV.X(); 
351     DComp1= DUV.X(); 
352     Up=Up1; 
353     Down=Down1;
354
355     while ( ( Abs(BSplS->UKnot(i)-Comp1)>Tol )&&(i!=Up1 ) )
356     {
357       i++;
358     }
359
360     cur = BSplS->UKnot(i); 
361   }
362   else if(Index==2)
363   {
364     i = Down2;
365     Comp1 = UV.Y();
366     DComp1=DUV.Y();
367     Up=Up2;
368     Down=Down2;
369
370     while ( ( Abs(BSplS->VKnot(i)-Comp1)>Tol )&&(i!=Up2 ) )
371     {
372       i++;
373     }
374
375     cur = BSplS->VKnot(i);
376   }
377   
378   if( Abs(Comp1-cur)<=Tol ) 
379   { 
380     Standard_Integer Bnd1 = Down, Bnd2 = Up;
381     if(Index==1)
382     {
383       TColStd_Array1OfReal Arr1(1,BSplS->NbUKnots());
384       BSplS->UKnots(Arr1); //   Up1=Arr1.Upper(); Down1=Arr1.Lower();
385       FindBounds(Arr1,cur,DUV.X(),Bnd1,Bnd2,DIsNull); 
386     }
387     else if(Index==2)
388     {
389       TColStd_Array1OfReal Arr2(1,BSplS->NbVKnots());
390       BSplS->VKnots(Arr2); //   Up2=Arr2.Upper(); Down2=Arr2.Lower();
391       FindBounds(Arr2,cur,DUV.Y(),Bnd1,Bnd2,DIsNull); 
392     }
393     
394     ReverseParam(Bnd1,Bnd2,Bnd1,Bnd2);
395
396     if(DIsNull==Standard_False)
397     {
398       if(Index==1)
399       {
400         LeftBot.SetX(BSplS->UKnot(Bnd1));
401         RightTop.SetX(BSplS->UKnot(Bnd2));
402       }
403       else if(Index==2)
404       { 
405         LeftBot.SetY(BSplS->VKnot(Bnd1));
406         RightTop.SetY(BSplS->VKnot(Bnd2));
407       }
408     }      
409   }
410   else//*********if Coord != Knot
411   {
412     if( (Index==1)&&(Comp1 < BSplS->UKnot(Down)) )
413     { 
414       LeftBot.SetX(BSplS->UKnot(Down)); 
415       RightTop.SetX( BSplS->UKnot(Down + 1) );
416       return;
417     }
418     else if( (Index==2)&&(Comp1 < BSplS->VKnot(Down)) )
419     { 
420       LeftBot.SetY(BSplS->VKnot(Down)); 
421       RightTop.SetY( BSplS->VKnot(Down + 1) );
422       return;
423     }
424     else if( (Index==1)&&(Comp1 > BSplS->UKnot(Up)) )
425     { 
426       RightTop.SetX(BSplS->UKnot(Up - 1)); 
427       LeftBot.SetX( BSplS->UKnot(Up) );
428       return;
429     }
430     else if( (Index==2)&&(Comp1 > BSplS->VKnot(Up)) )
431     { 
432       RightTop.SetY(BSplS->VKnot(Up - 1)); 
433       LeftBot.SetY( BSplS->VKnot(Up) );
434       return;
435     } 
436     else
437     {
438       Standard_Real f = 0., l = 1.;
439       if (Index==1)
440       {
441         f=BSplS->UKnot(Down);
442         l=BSplS->UKnot(Up);
443       }
444       else if (Index == 2)
445       {
446         f=BSplS->VKnot(Down);
447         l=BSplS->VKnot(Up);
448       }
449
450       i = Down;
451       if ((!(Comp1 < f))&&(!(Comp1 > l)))
452       {
453         if (Index==1)
454         {
455           while (!(((f=BSplS->UKnot(i)) < Comp1)&&((l=BSplS->UKnot(i+1)) > Comp1)) && (i<Up)) 
456           {
457             i++;
458           }
459         }
460         else if (Index==2)
461         {
462           while (!(((f=BSplS->VKnot(i)) < Comp1)&&((l=BSplS->VKnot(i+1)) > Comp1)) && (i<Up)) 
463           {
464             i++;
465           }
466         }
467       }
468       else
469         ReverseParam(f,l,f,l);
470
471       if(i!=Up)
472       {
473         if(Abs(DComp1)>Tol)
474         {
475           if(Index==1) 
476           {  
477             if(DComp1>0)
478             { 
479               LeftBot.SetX(Comp1); 
480               RightTop.SetX(l);
481             }
482             else if(DComp1<0)
483             {
484               LeftBot.SetX(f); 
485               RightTop.SetX(Comp1);
486             }
487           }
488           else if(Index==2)
489           {
490             if(DComp1>0)
491             { 
492               LeftBot.SetY(Comp1); 
493               RightTop.SetY(l);
494             }
495             else if(DComp1<0)
496             { 
497               LeftBot.SetY(f); 
498               RightTop.SetY(Comp1);
499             }
500           }
501         }
502         else
503         {
504           if(Abs(DComp1)<Tol)
505           {
506             if(Index==1)
507             {
508               LeftBot.SetX(f); 
509               RightTop.SetX(l);
510             }
511             else if(Index==2) 
512             { 
513               LeftBot.SetY(f); 
514               RightTop.SetY(l);
515             }
516           }
517         }
518       } 
519       else if(i==Up)
520       {
521         if(Index==1)
522         {
523           LeftBot.SetX(Comp1);
524           RightTop.SetX(BSplS->UKnot(i));
525         }
526         else if(Index==2)
527         { 
528           LeftBot.SetY(Comp1); 
529           RightTop.SetY(BSplS->VKnot(i));
530         }
531       }
532     }
533   }
534 }
535 //=======================================================================
536 //function :Locate2Coord
537 //purpose  : along non-BSpline curve
538 //=======================================================================
539
540
541 static void Locate2Coord(const Standard_Integer Index,
542                          const gp_Pnt2d& UV, const gp_Vec2d& DUV, 
543                          const Standard_Real I1, 
544                          const Standard_Real I2,
545                          gp_Pnt2d& LeftBot, gp_Pnt2d& RightTop)
546 {
547   Standard_Real Tol=Precision::PConfusion()/10;
548   Standard_Real Comp1=0,DComp1=0;
549   if(Index==1)   {   Comp1=UV.X();
550                      DComp1=DUV.X();} 
551   else
552     if(Index==2)   {Comp1=UV.Y();
553                     DComp1=DUV.Y();}
554   
555   if((Comp1!=I1)&&(Comp1!=I2))
556     { if(Abs(DComp1) > Tol)
557         {  if(DComp1 < 0) 
558              {  if(Index==1) {  LeftBot.SetX(I1);
559                                 RightTop.SetX(Comp1);}
560                 if(Index==2) {  LeftBot.SetY(I1);
561                                 RightTop.SetY(Comp1);}
562               }
563         else
564           if(DComp1 > 0) 
565             {   if(Index==1) {  LeftBot.SetX(Comp1);
566                                 RightTop.SetX(I2);}
567                 if(Index==2) {  LeftBot.SetY(Comp1);
568                                 RightTop.SetY(I2);}
569               }
570           else { if(Index==1) { LeftBot.SetX(I1); 
571                                 RightTop.SetX(I2);}
572                  if(Index==2) { LeftBot.SetY(I1);
573                                 RightTop.SetY(I2);}
574                }
575          }
576     else 
577       if(Abs(DComp1)<=Tol) {
578                               if(Index==1) { LeftBot.SetX(I1) ;
579                                              RightTop.SetX(I2);}
580                               if(Index==2) { LeftBot.SetY(I1) ;
581                                              RightTop.SetY(I2);}
582                             }
583     }else 
584       if(Abs(Comp1-I1)<Tol)
585         { if(Index==1) { LeftBot.SetX(I1) ;
586                          RightTop.SetX(I2);}
587           if(Index==2) { LeftBot.SetY(I1) ;
588                          RightTop.SetY(I2);}
589         }
590       else
591         if(Abs(Comp1-I2)<Tol)
592           {        if(Index==1) { LeftBot.SetX(I1); 
593                                   RightTop.SetX(I2);}
594                    if(Index==2) { LeftBot.SetY(I1);
595                                   RightTop.SetY(I2);}
596                  } 
597 }
598
599 //=======================================================================
600 //function :Locate2Coord
601 //purpose  : 
602 //=======================================================================
603
604 static void Locate2Coord(const Standard_Integer Index,
605                          const gp_Pnt2d& UV, const gp_Vec2d& DUV, 
606                          const Handle(Geom_BSplineSurface)& BSplS,
607                          const TColStd_Array1OfReal& Arr,
608                          gp_Pnt2d& LeftBot, gp_Pnt2d& RightTop)
609 {
610   Standard_Real Comp=0,DComp=0,Tmp1=0.0,Tmp2=0.0;
611   Standard_Real Tol=Precision::PConfusion()/10;
612   Standard_Integer N=0, NUp=0, NLo=0;
613   if(Index==1)
614     { Comp=UV.X();
615       DComp=DUV.Y();
616       NUp = BSplS->LastUKnotIndex();
617       NLo = BSplS->FirstUKnotIndex();
618     }
619   if(Index==2)
620      { Comp=UV.Y();
621        DComp=DUV.X(); 
622        NUp = BSplS->LastVKnotIndex();
623        NLo = BSplS->FirstVKnotIndex();
624      }
625   
626   if((DComp > 0)&&(Abs(DComp)>Tol)) {  
627     Hunt(Arr,Comp,N);
628     if (N >= NUp){
629       //limit case: Hunt() cought upper knot. Take the last span. 
630       N = NUp - 1;
631     }
632     if(Index==1) {  Tmp1=BSplS->UKnot(N);
633                     Tmp2=BSplS->UKnot(N+1);}
634     else
635       if(Index==2) {      Tmp1=BSplS->VKnot(N);
636                           Tmp2=BSplS->VKnot(N+1);}
637     
638     ReverseParam(Tmp1,Tmp2,Tmp1,Tmp2);
639     
640     if(Index==1) {      LeftBot.SetX(Tmp1);
641                         RightTop.SetX(Tmp2);}
642     else
643       if(Index==2) {      LeftBot.SetY(Tmp1);
644                           RightTop.SetY(Tmp2);}
645     }
646   else
647     if((DComp < 0)&&(Abs(DComp)>Tol)){      
648      Hunt(Arr,Comp,N);
649      if (N <= NLo) {
650        //limit case: Hunt() cought lower knot. Take the first span.
651        N = NLo + 1;
652      }
653      if(Index==1) {       Tmp1=BSplS->UKnot(N-1);
654                           Tmp2=BSplS->UKnot(N);}
655         else
656           if(Index==2) { Tmp1=BSplS->VKnot(N-1);
657                          Tmp2=BSplS->VKnot(N);}
658      
659         ReverseParam(Tmp1,Tmp2,Tmp1,Tmp2);
660      
661         if(Index==1) {    LeftBot.SetX(Tmp1);
662                           RightTop.SetX(Tmp2);}
663         else
664           if(Index==2) {            LeftBot.SetY(Tmp1);
665                                     RightTop.SetY(Tmp2);}
666    } 
667 }
668
669
670
671 //=======================================================================
672 //function : Adaptor3d_CurveOnSurface
673 //purpose  : 
674 //=======================================================================
675
676 Adaptor3d_CurveOnSurface::Adaptor3d_CurveOnSurface()
677      : myType(GeomAbs_OtherCurve), myIntCont(GeomAbs_CN)
678 {}
679
680 //=======================================================================
681 //function : Adaptor3d_CurveOnSurface
682 //purpose  : 
683 //=======================================================================
684
685 Adaptor3d_CurveOnSurface::Adaptor3d_CurveOnSurface
686 (const Handle(Adaptor3d_HSurface)& S) 
687      : myType(GeomAbs_OtherCurve), myIntCont(GeomAbs_CN)
688 {
689   Load(S);
690 }
691
692 //=======================================================================
693 //function : Adaptor3d_CurveOnSurface
694 //purpose  : 
695 //=======================================================================
696
697 Adaptor3d_CurveOnSurface::Adaptor3d_CurveOnSurface
698 (const Handle(Adaptor2d_HCurve2d)& C,
699  const Handle(Adaptor3d_HSurface)& S)
700      : myType(GeomAbs_OtherCurve), myIntCont(GeomAbs_CN)
701 {
702   Load(S);
703   Load(C);
704 }
705
706 //=======================================================================
707 //function : Load
708 //purpose  : 
709 //=======================================================================
710
711 void Adaptor3d_CurveOnSurface::Load(const Handle(Adaptor3d_HSurface)& S) 
712 {
713   mySurface = S;
714   if (!myCurve.IsNull()) EvalKPart();
715 }
716
717 //=======================================================================
718 //function : Load
719 //purpose  : 
720 //=======================================================================
721
722 void Adaptor3d_CurveOnSurface::Load(const Handle(Adaptor2d_HCurve2d)& C) 
723 {
724   myCurve = C;
725   if (!mySurface.IsNull()) 
726      {
727        EvalKPart();
728        GeomAbs_SurfaceType  SType ;
729        SType = mySurface->GetType();
730        if( SType == GeomAbs_BSplineSurface)
731          EvalFirstLastSurf();
732        if( SType == GeomAbs_SurfaceOfExtrusion)
733          EvalFirstLastSurf();
734          if( SType == GeomAbs_SurfaceOfRevolution)
735          EvalFirstLastSurf();
736          if( SType == GeomAbs_OffsetSurface) {
737           SType = mySurface->BasisSurface()->GetType();
738           if( SType == GeomAbs_SurfaceOfRevolution ||
739               SType == GeomAbs_SurfaceOfExtrusion ||
740               SType == GeomAbs_BSplineSurface )
741             EvalFirstLastSurf();
742         }
743      }
744 }
745
746 //=======================================================================
747 //function : FirstParameter
748 //purpose  : 
749 //=======================================================================
750
751 Standard_Real Adaptor3d_CurveOnSurface::FirstParameter() const 
752 {
753   return myCurve->FirstParameter();
754 }
755
756 //=======================================================================
757 //function : LastParameter
758 //purpose  : 
759 //=======================================================================
760
761 Standard_Real Adaptor3d_CurveOnSurface::LastParameter() const 
762 {
763   return myCurve->LastParameter();
764 }
765
766 //=======================================================================
767 //function : Continuity
768 //purpose  : 
769 //=======================================================================
770
771 GeomAbs_Shape Adaptor3d_CurveOnSurface::Continuity() const
772 {
773   GeomAbs_Shape ContC  = myCurve->Continuity();
774   GeomAbs_Shape ContSu = mySurface->UContinuity();
775   if ( ContSu < ContC) ContC = ContSu;
776   GeomAbs_Shape ContSv = mySurface->VContinuity();
777   if ( ContSv < ContC) ContC = ContSv;
778
779   return ContC;
780 }
781
782 // Auxiliary: adds roots of equation to sorted sequence of parameters
783 // along curve, keeping it sorted and avoiding repetitions (within tolerance Tol)
784 static void AddIntervals (const Handle(TColStd_HSequenceOfReal)& theParameters, 
785                           const math_FunctionRoots& theRoots, Standard_Real theTol)
786 {
787   if (! theRoots.IsDone() || theRoots.IsAllNull())
788     return;
789
790   Standard_Integer nsol = theRoots.NbSolutions();
791   for (Standard_Integer i = 1; i <= nsol; i++)
792   {
793     Standard_Real param = theRoots.Value(i);
794     if (param - theParameters->Value(1) < theTol) // skip param if equal to or less than theParameters(1)
795       continue;
796     for (Standard_Integer j=2; j <= theParameters->Length(); ++j)
797     {
798       Standard_Real aDelta = theParameters->Value(j) - param;
799       if (aDelta > theTol)
800       {
801         theParameters->InsertBefore (j, param);
802         break;
803       }
804       else if (aDelta >= -theTol) // param == theParameters(j) within Tol
805         break;
806     }
807   }
808 }
809
810 //=======================================================================
811 //function : NbIntervals
812 //purpose  : 
813 //=======================================================================
814
815 Standard_Integer Adaptor3d_CurveOnSurface::NbIntervals (const GeomAbs_Shape S) const
816 {
817   if(S == myIntCont && !myIntervals.IsNull())
818     return myIntervals->Length()-1;
819   
820   Standard_Integer nu,nv,nc;
821   nu=mySurface->NbUIntervals(S);
822   nv=mySurface->NbVIntervals(S);
823
824   TColStd_Array1OfReal TabU(1,nu+1);
825   TColStd_Array1OfReal TabV(1,nv+1);
826   Standard_Integer NbSample = 20;
827   Standard_Real U,V,Tdeb,Tfin;
828   Tdeb=myCurve->FirstParameter();
829   Tfin=myCurve->LastParameter();
830
831   nc=myCurve->NbIntervals(S);
832   TColStd_Array1OfReal TabC(1,nc+1);
833   myCurve->Intervals(TabC,S);
834
835   Standard_Real Tol= Precision::PConfusion()/10;
836
837   // sorted sequence of parameters defining continuity intervals;
838   // started with own intervals of curve and completed by 
839   // additional points coming from surface discontinuities
840   Handle(TColStd_HSequenceOfReal) aIntervals = new TColStd_HSequenceOfReal;
841   for (Standard_Integer i = 1; i <= nc + 1; i++)
842   {
843     aIntervals->Append(TabC(i));
844   }
845  
846   if (nu>1)
847   {
848     mySurface->UIntervals(TabU,S);
849     for(Standard_Integer iu = 2;iu <= nu; iu++)
850     {
851       U = TabU.Value(iu);
852       Adaptor3d_InterFunc Func(myCurve,U,1);
853       math_FunctionRoots Resol(Func,Tdeb,Tfin,NbSample,Tol,Tol,Tol,0.);
854       AddIntervals (aIntervals, Resol, Tol);
855     }
856   }
857   if (nv>1)
858   {
859     mySurface->VIntervals(TabV,S);
860     for(Standard_Integer iv = 2;iv <= nv; iv++)
861     {
862       V = TabV.Value(iv);
863       Adaptor3d_InterFunc Func(myCurve,V,2);
864       math_FunctionRoots Resol(Func,Tdeb,Tfin,NbSample,Tol,Tol,Tol,0.);
865       AddIntervals (aIntervals, Resol, Tol);
866     }
867   }
868
869   // for case intervals==1 and first point == last point SequenceOfReal
870   // contains only one value, therefore it is necessary to add second
871   // value into aIntervals which will be equal first value.
872   if (aIntervals->Length() == 1)
873     aIntervals->Append (aIntervals->Value(1));
874
875   const_cast<Adaptor3d_CurveOnSurface*>(this)->myIntervals = aIntervals;
876   const_cast<Adaptor3d_CurveOnSurface*>(this)->myIntCont = S;
877   return myIntervals->Length() - 1;
878 }
879
880 //=======================================================================
881 //function : Intervals
882 //purpose  : 
883 //=======================================================================
884
885 void Adaptor3d_CurveOnSurface::Intervals(TColStd_Array1OfReal& T,
886                                          const GeomAbs_Shape S) const
887 {
888   NbIntervals(S);
889   Standard_ASSERT_RAISE (T.Length() == myIntervals->Length(), "Error: Wrong size of array buffer in call to Adaptor3d_CurveOnSurface::Intervals");
890   for(Standard_Integer i=1; i<=myIntervals->Length(); i++) {
891     T(i) = myIntervals->Value(i);
892   }
893 }
894
895 //=======================================================================
896 //function : Trim
897 //purpose  : 
898 //=======================================================================
899
900 Handle(Adaptor3d_HCurve) Adaptor3d_CurveOnSurface::Trim
901 (const Standard_Real First, 
902  const Standard_Real Last, 
903  const Standard_Real Tol) const 
904 {
905   Handle(Adaptor3d_HCurveOnSurface) HCS = new Adaptor3d_HCurveOnSurface();
906   HCS->ChangeCurve().Load(mySurface);
907   HCS->ChangeCurve().Load(myCurve->Trim(First,Last,Tol));
908   return HCS;
909 }
910
911 //=======================================================================
912 //function : IsClosed
913 //purpose  : 
914 //=======================================================================
915
916 Standard_Boolean Adaptor3d_CurveOnSurface::IsClosed() const
917 {
918   return myCurve->IsClosed();
919 }
920
921 //=======================================================================
922 //function : IsPeriodic
923 //purpose  : 
924 //=======================================================================
925
926 Standard_Boolean Adaptor3d_CurveOnSurface::IsPeriodic() const
927 {
928   return myCurve->IsPeriodic();
929 }
930
931 //=======================================================================
932 //function : Period
933 //purpose  : 
934 //=======================================================================
935
936 Standard_Real Adaptor3d_CurveOnSurface::Period() const
937 {
938   return myCurve->Period();
939 }
940
941 //=======================================================================
942 //function : Value
943 //purpose  : 
944 //=======================================================================
945
946 gp_Pnt Adaptor3d_CurveOnSurface::Value(const Standard_Real U ) const
947 {
948   gp_Pnt P;
949   gp_Pnt2d Puv;
950   
951   if      (myType == GeomAbs_Line  ) P = ElCLib::Value(U,myLin );
952   else if (myType == GeomAbs_Circle) P = ElCLib::Value(U,myCirc);
953   else {
954     myCurve->D0(U,Puv);
955     mySurface->D0(Puv.X(),Puv.Y(),P);
956   }
957
958   return P;
959 }
960
961 //=======================================================================
962 //function : D0
963 //purpose  : 
964 //=======================================================================
965
966 void Adaptor3d_CurveOnSurface::D0(const Standard_Real U ,
967                                   gp_Pnt& P) const 
968 {
969   gp_Pnt2d Puv;
970   
971   if      (myType == GeomAbs_Line  ) P = ElCLib::Value(U,myLin );
972   else if (myType == GeomAbs_Circle) P = ElCLib::Value(U,myCirc);
973   else {
974     myCurve->D0(U,Puv);
975     mySurface->D0(Puv.X(),Puv.Y(),P);
976   }
977
978 }
979
980
981 //=======================================================================
982 //function : D1
983 //purpose  : 
984 //=======================================================================
985
986 void Adaptor3d_CurveOnSurface::D1(const Standard_Real U ,
987                                 gp_Pnt& P,
988                                 gp_Vec& V) const 
989 {
990   gp_Pnt2d Puv;
991   gp_Vec2d Duv;
992   gp_Vec D1U,D1V;
993   
994   Standard_Real FP = myCurve->FirstParameter();
995   Standard_Real LP = myCurve->LastParameter();
996   
997   Standard_Real Tol= Precision::PConfusion()/10; 
998   if( ( Abs(U-FP)<Tol)&&(!myFirstSurf.IsNull()) )
999     {
1000       myCurve->D1(U,Puv,Duv);
1001       myFirstSurf->D1(Puv.X(),Puv.Y(),P,D1U,D1V);
1002       V.SetLinearForm(Duv.X(),D1U,Duv.Y(),D1V);
1003     }
1004   else
1005     if( (Abs(U-LP)<Tol)&&(!myLastSurf.IsNull()) )
1006       {
1007         myCurve->D1(U,Puv,Duv);
1008         myLastSurf->D1(Puv.X(),Puv.Y(),P,D1U,D1V);
1009         V.SetLinearForm(Duv.X(),D1U,Duv.Y(),D1V);
1010       }
1011     else
1012       if      (myType == GeomAbs_Line  ) ElCLib::D1(U,myLin ,P,V);
1013       else if (myType == GeomAbs_Circle) ElCLib::D1(U,myCirc,P,V);
1014       else {
1015         myCurve->D1(U,Puv,Duv);
1016         mySurface->D1(Puv.X(),Puv.Y(),P,D1U,D1V);
1017         V.SetLinearForm(Duv.X(),D1U,Duv.Y(),D1V);
1018       }
1019 }
1020 //=======================================================================
1021 //function : D2
1022 //purpose  : 
1023 //=======================================================================
1024
1025 void Adaptor3d_CurveOnSurface::D2(const Standard_Real U,
1026                                 gp_Pnt& P,
1027                                 gp_Vec& V1,
1028                                 gp_Vec& V2) const
1029
1030   gp_Pnt2d UV;
1031   gp_Vec2d DW,D2W; 
1032   gp_Vec D1U,D1V,D2U,D2V,D2UV;
1033   
1034   Standard_Real FP = myCurve->FirstParameter();
1035   Standard_Real LP = myCurve->LastParameter();
1036   
1037   Standard_Real Tol= Precision::PConfusion()/10; 
1038   if( (Abs(U-FP)<Tol)&&(!myFirstSurf.IsNull()) )
1039     {
1040       myCurve->D2(U,UV,DW,D2W);
1041       myFirstSurf->D2(UV.X(),UV.Y(),P,D1U,D1V,D2U,D2V,D2UV);
1042       
1043       V1.SetLinearForm(DW.X(),D1U,DW.Y(),D1V);
1044       V2.SetLinearForm(D2W.X(), D1U, D2W.Y(), D1V, 2.*DW.X()*DW.Y(),D2UV);
1045       V2.SetLinearForm(DW.X()*DW.X(), D2U, DW.Y()*DW.Y(), D2V, V2);  
1046     }
1047   else  
1048     if( (Abs(U-LP)<Tol)&&(!myLastSurf.IsNull()) )
1049       {
1050         myCurve->D2(U,UV,DW,D2W);
1051         myLastSurf->D2(UV.X(),UV.Y(),P,D1U,D1V,D2U,D2V,D2UV);
1052         
1053         V1.SetLinearForm(DW.X(),D1U,DW.Y(),D1V);
1054         V2.SetLinearForm(D2W.X(), D1U, D2W.Y(), D1V, 2.*DW.X()*DW.Y(),D2UV);
1055         V2.SetLinearForm(DW.X()*DW.X(), D2U, DW.Y()*DW.Y(), D2V, V2);  
1056       }
1057     else
1058       if (myType == GeomAbs_Line  ) {
1059         ElCLib::D1(U,myLin,P,V1);
1060           V2.SetCoord(0.,0.,0.);
1061         }
1062       else if (myType == GeomAbs_Circle) ElCLib::D2(U,myCirc,P,V1,V2);
1063       else {
1064         myCurve->D2(U,UV,DW,D2W);
1065         mySurface->D2(UV.X(),UV.Y(),P,D1U,D1V,D2U,D2V,D2UV);
1066         
1067         V1.SetLinearForm(DW.X(),D1U,DW.Y(),D1V);
1068         V2.SetLinearForm(D2W.X(), D1U, D2W.Y(), D1V, 2.*DW.X()*DW.Y(),D2UV);
1069         V2.SetLinearForm(DW.X()*DW.X(), D2U, DW.Y()*DW.Y(), D2V, V2);  
1070       }
1071 }
1072
1073 //=======================================================================
1074 //function : D3
1075 //purpose  : 
1076 //=======================================================================
1077
1078 void Adaptor3d_CurveOnSurface::D3
1079   (const Standard_Real U,
1080    gp_Pnt& P,
1081    gp_Vec& V1,
1082    gp_Vec& V2,
1083    gp_Vec& V3) const
1084
1085   
1086   Standard_Real Tol= Precision::PConfusion()/10; 
1087   gp_Pnt2d UV;
1088   gp_Vec2d DW,D2W,D3W;
1089   gp_Vec D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV;
1090   
1091   Standard_Real FP = myCurve->FirstParameter();
1092   Standard_Real LP = myCurve->LastParameter();
1093   
1094   if( (Abs(U-FP)<Tol)&&(!myFirstSurf.IsNull()) )
1095     { myCurve->D3(U,UV,DW,D2W,D3W);
1096       myFirstSurf->D3(UV.X(),UV.Y(),P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
1097       V1.SetLinearForm(DW.X(),D1U,DW.Y(),D1V);
1098       V2.SetLinearForm(D2W.X(), D1U, D2W.Y(), D1V, 2.*DW.X()*DW.Y(),D2UV);
1099       V2.SetLinearForm(DW.X()*DW.X(), D2U, DW.Y()*DW.Y(), D2V, V2); 
1100       V3=SetLinearForm( DW, D2W, D3W, D1U, D1V, D2U, D2V, D2UV, D3U, D3V, D3UUV, D3UVV);
1101     }else
1102       
1103       if( (Abs(U-LP)<Tol)&&(!myLastSurf.IsNull()) )
1104         { myCurve->D3(U,UV,DW,D2W,D3W);
1105           myLastSurf->D3(UV.X(),UV.Y(),P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
1106           V1.SetLinearForm(DW.X(),D1U,DW.Y(),D1V);
1107           
1108           V2.SetLinearForm(D2W.X(), D1U, D2W.Y(), D1V, 2.*DW.X()*DW.Y(),D2UV);
1109           V2.SetLinearForm(DW.X()*DW.X(), D2U, DW.Y()*DW.Y(), D2V, V2);  
1110           V3=SetLinearForm( DW, D2W, D3W, D1U, D1V, D2U, D2V, D2UV, D3U, D3V, D3UUV, D3UVV);  
1111         }
1112       else
1113         if (myType == GeomAbs_Line  ) {
1114           ElCLib::D1(U,myLin,P,V1);
1115             V2.SetCoord(0.,0.,0.);
1116             V3.SetCoord(0.,0.,0.);
1117           }
1118         else if (myType == GeomAbs_Circle) ElCLib::D3(U,myCirc,P,V1,V2,V3);
1119         else {
1120           myCurve->D3(U,UV,DW,D2W,D3W);
1121           mySurface->D3(UV.X(),UV.Y(),P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
1122           V1.SetLinearForm(DW.X(),D1U,DW.Y(),D1V);
1123           
1124           V2.SetLinearForm(D2W.X(), D1U, D2W.Y(), D1V, 2.*DW.X()*DW.Y(),D2UV);
1125           V2.SetLinearForm(DW.X()*DW.X(), D2U, DW.Y()*DW.Y(), D2V, V2);  
1126           V3=SetLinearForm( DW, D2W, D3W, D1U, D1V, D2U, D2V, D2UV, D3U, D3V, D3UUV, D3UVV);  
1127         }
1128 }
1129
1130
1131 //=======================================================================
1132 //function : DN
1133 //purpose  : 
1134 //=======================================================================
1135
1136 gp_Vec Adaptor3d_CurveOnSurface::DN
1137   (const Standard_Real U, 
1138    const Standard_Integer N) const 
1139 {
1140   gp_Pnt P;
1141   gp_Vec V1, V2, V;
1142   switch (N) {
1143   case 1:
1144     D1(U,P,V);
1145     break ;
1146   case 2: 
1147     D2(U,P,V1,V);
1148     break ;
1149   case 3: 
1150     D3(U,P,V1,V2,V);
1151     break ;
1152   default:
1153     Standard_NotImplemented::Raise("Adaptor3d_CurveOnSurface:DN");
1154     break;
1155   }
1156   return V;
1157 }
1158
1159
1160 //=======================================================================
1161 //function : Resolution
1162 //purpose  : 
1163 //=======================================================================
1164
1165 Standard_Real Adaptor3d_CurveOnSurface::Resolution
1166   (const Standard_Real R3d) const
1167 {
1168   Standard_Real ru,rv;
1169   ru = mySurface->UResolution(R3d);
1170   rv = mySurface->VResolution(R3d);
1171   return myCurve->Resolution(Min(ru,rv)); 
1172 }
1173
1174
1175 //=======================================================================
1176 //function : GetType
1177 //purpose  : 
1178 //=======================================================================
1179
1180 GeomAbs_CurveType Adaptor3d_CurveOnSurface::GetType() const 
1181 {
1182   return myType;
1183
1184
1185
1186 //=======================================================================
1187 //function : Line
1188 //purpose  : 
1189 //=======================================================================
1190
1191 gp_Lin Adaptor3d_CurveOnSurface::Line() const
1192 {
1193   Standard_NoSuchObject_Raise_if(myType != GeomAbs_Line, "Adaptor3d_CurveOnSurface::Line(): curve is not a line")
1194   return myLin;
1195 }
1196
1197 //=======================================================================
1198 //function : Circle
1199 //purpose  : 
1200 //=======================================================================
1201
1202 gp_Circ Adaptor3d_CurveOnSurface::Circle() const
1203 {
1204   Standard_NoSuchObject_Raise_if(myType != GeomAbs_Circle, "Adaptor3d_CurveOnSurface::Line(): curve is not a circle")
1205   return myCirc;
1206 }
1207
1208 //=======================================================================
1209 //function : Ellipse
1210 //purpose  : 
1211 //=======================================================================
1212
1213 gp_Elips Adaptor3d_CurveOnSurface::Ellipse() const
1214 {
1215   return to3d(mySurface->Plane(),myCurve->Ellipse());
1216 }
1217
1218 //=======================================================================
1219 //function : Hyperbola
1220 //purpose  : 
1221 //=======================================================================
1222
1223 gp_Hypr Adaptor3d_CurveOnSurface::Hyperbola() const
1224 {
1225   return to3d(mySurface->Plane(),myCurve->Hyperbola());
1226 }
1227
1228 //=======================================================================
1229 //function : Parabola
1230 //purpose  : 
1231 //=======================================================================
1232
1233 gp_Parab Adaptor3d_CurveOnSurface::Parabola() const
1234 {
1235   return to3d(mySurface->Plane(),myCurve->Parabola());
1236 }
1237
1238 Standard_Integer  Adaptor3d_CurveOnSurface::Degree() const
1239 {
1240   
1241   // on a parametric surface should multiply
1242   // return TheCurve2dTool::Degree(myCurve);
1243
1244   return myCurve->Degree();
1245 }
1246
1247 //=======================================================================
1248 //function : IsRational
1249 //purpose  : 
1250 //=======================================================================
1251
1252 Standard_Boolean Adaptor3d_CurveOnSurface::IsRational() const
1253 {
1254   return ( myCurve->IsRational() ||
1255           mySurface->IsURational() ||
1256           mySurface->IsVRational()   );
1257 }
1258
1259 //=======================================================================
1260 //function : NbPoles
1261 //purpose  : 
1262 //=======================================================================
1263
1264 Standard_Integer  Adaptor3d_CurveOnSurface::NbPoles() const
1265 {
1266   // on a parametric surface should multiply
1267   return myCurve->NbPoles();
1268 }
1269
1270 //=======================================================================
1271 //function : NbKnots
1272 //purpose  : 
1273 //=======================================================================
1274
1275 Standard_Integer Adaptor3d_CurveOnSurface::NbKnots() const {
1276   if (mySurface->GetType()==GeomAbs_Plane)
1277     return myCurve->NbKnots();
1278   else {
1279     Standard_NoSuchObject::Raise();
1280     return 0;
1281   }
1282 }
1283
1284 //=======================================================================
1285 //function : Bezier
1286 //purpose  : 
1287 //=======================================================================
1288
1289 Handle(Geom_BezierCurve) Adaptor3d_CurveOnSurface::Bezier() const  
1290 {
1291   Standard_NoSuchObject_Raise_if
1292     ( mySurface->GetType() != GeomAbs_Plane,
1293      "Adaptor3d_CurveOnSurface : Bezier");
1294
1295   Handle(Geom2d_BezierCurve) Bez2d = myCurve->Bezier();
1296   Standard_Integer NbPoles = Bez2d->NbPoles();
1297
1298   const gp_Pln& Plane = mySurface->Plane();
1299
1300   TColgp_Array1OfPnt Poles(1,NbPoles);
1301   for ( Standard_Integer i=1; i<= NbPoles; i++) {
1302     Poles(i) = to3d( Plane, Bez2d->Pole(i));
1303   }
1304   Handle(Geom_BezierCurve) Bez;
1305
1306   if (Bez2d->IsRational()) {
1307     TColStd_Array1OfReal Weights(1,NbPoles);
1308     Bez2d->Weights(Weights);
1309     Bez = new Geom_BezierCurve(Poles,Weights);
1310   }
1311   else {
1312     Bez = new Geom_BezierCurve(Poles);
1313   }
1314   return Bez;
1315 }
1316
1317 //=======================================================================
1318 //function : BSpline
1319 //purpose  : 
1320 //=======================================================================
1321
1322 Handle(Geom_BSplineCurve) Adaptor3d_CurveOnSurface::BSpline() const  
1323 {
1324   Standard_NoSuchObject_Raise_if
1325     ( mySurface->GetType() != GeomAbs_Plane,
1326      "Adaptor3d_CurveOnSurface : BSpline");
1327
1328   Handle(Geom2d_BSplineCurve) Bsp2d = myCurve->BSpline();
1329   Standard_Integer NbPoles = Bsp2d->NbPoles();
1330
1331   const gp_Pln& Plane = mySurface->Plane();
1332
1333   TColgp_Array1OfPnt Poles(1,NbPoles);
1334   for ( Standard_Integer i=1; i<= NbPoles; i++) {
1335     Poles(i) = to3d( Plane, Bsp2d->Pole(i));
1336   }
1337
1338   TColStd_Array1OfReal    Knots(1,Bsp2d->NbKnots());
1339   TColStd_Array1OfInteger Mults(1,Bsp2d->NbKnots());
1340   Bsp2d->Knots(Knots);
1341   Bsp2d->Multiplicities(Mults);
1342
1343   Handle(Geom_BSplineCurve) Bsp;
1344
1345   if (Bsp2d->IsRational()) {
1346     TColStd_Array1OfReal Weights(1,NbPoles);
1347     Bsp2d->Weights(Weights);
1348     Bsp = new Geom_BSplineCurve(Poles,Weights,Knots,Mults,
1349                                 Bsp2d->Degree(),
1350                                 Bsp2d->IsPeriodic());
1351   }
1352   else {
1353     Bsp = new Geom_BSplineCurve(Poles,Knots,Mults,
1354                                 Bsp2d->Degree(),
1355                                 Bsp2d->IsPeriodic());
1356   }
1357   return Bsp;
1358 }
1359
1360 //=======================================================================
1361 //function : GetCurve
1362 //purpose  : 
1363 //=======================================================================
1364
1365 const Handle(Adaptor2d_HCurve2d)& Adaptor3d_CurveOnSurface::GetCurve() const
1366 {
1367   return myCurve;
1368 }
1369
1370 //=======================================================================
1371 //function : GetSurface
1372 //purpose  : 
1373 //=======================================================================
1374
1375 const Handle(Adaptor3d_HSurface)& Adaptor3d_CurveOnSurface::GetSurface() const 
1376 {
1377   return mySurface;
1378 }
1379
1380 //=======================================================================
1381 //function : ChangeCurve
1382 //purpose  : 
1383 //=======================================================================
1384
1385 Handle(Adaptor2d_HCurve2d)& Adaptor3d_CurveOnSurface::ChangeCurve() 
1386 {
1387   return myCurve;
1388 }
1389
1390 //=======================================================================
1391 //function : ChangeSurface
1392 //purpose  : 
1393 //=======================================================================
1394
1395 Handle(Adaptor3d_HSurface)& Adaptor3d_CurveOnSurface::ChangeSurface() {
1396   return mySurface;
1397 }
1398
1399 //=======================================================================
1400 //function : EvalKPart
1401 //purpose  : 
1402 //=======================================================================
1403
1404 void Adaptor3d_CurveOnSurface::EvalKPart()
1405 {
1406   myType = GeomAbs_OtherCurve;
1407
1408   GeomAbs_SurfaceType STy = mySurface->GetType();
1409   GeomAbs_CurveType   CTy = myCurve->GetType();
1410   if (STy == GeomAbs_Plane) {
1411     myType =  CTy;
1412     if (myType == GeomAbs_Circle) 
1413       myCirc = to3d(mySurface->Plane(),myCurve->Circle());
1414     else if (myType == GeomAbs_Line) {
1415       gp_Pnt P;
1416       gp_Vec V;
1417       gp_Pnt2d Puv;
1418       gp_Vec2d Duv;
1419       myCurve->D1(0.,Puv,Duv);
1420       gp_Vec D1U,D1V;
1421       mySurface->D1(Puv.X(),Puv.Y(),P,D1U,D1V);
1422       V.SetLinearForm(Duv.X(),D1U,Duv.Y(),D1V); 
1423       myLin = gp_Lin(P,V);
1424     }
1425   }
1426   else {
1427     if ( CTy == GeomAbs_Line) {
1428       gp_Dir2d D = myCurve->Line().Direction();
1429       if ( D.IsParallel(gp::DX2d(),Precision::Angular())) { // Iso V.
1430         if ( STy == GeomAbs_Sphere) {
1431           gp_Pnt2d  P    = myCurve->Line().Location();
1432           if ( Abs( Abs(P.Y()) -M_PI/2. ) >= Precision::PConfusion()) {
1433             myType = GeomAbs_Circle;
1434             gp_Sphere Sph  =  mySurface->Sphere();
1435             gp_Ax3    Axis = Sph.Position();
1436             myCirc   = ElSLib::SphereVIso(Axis,
1437                                           Sph.Radius(),
1438                                           P.Y());
1439             gp_Dir DRev = Axis.XDirection().Crossed(Axis.YDirection());
1440             gp_Ax1 AxeRev(Axis.Location(), DRev);
1441             myCirc.Rotate(AxeRev, P.X());
1442             if ( D.IsOpposite(gp::DX2d(),Precision::Angular())) {
1443               gp_Ax2 Ax = myCirc.Position();
1444               Ax.SetDirection(Ax.Direction().Reversed());
1445               myCirc.SetPosition(Ax);
1446             }
1447           }
1448         }
1449         else if ( STy == GeomAbs_Cylinder) {
1450           myType = GeomAbs_Circle;
1451           gp_Cylinder Cyl  = mySurface->Cylinder();
1452           gp_Pnt2d    P    = myCurve->Line().Location();
1453           gp_Ax3      Axis = Cyl.Position();
1454           myCirc           = ElSLib::CylinderVIso(Axis,
1455                                                   Cyl.Radius(),
1456                                                   P.Y());
1457           gp_Dir DRev = Axis.XDirection().Crossed(Axis.YDirection());
1458           gp_Ax1 AxeRev(Axis.Location(), DRev);
1459           myCirc.Rotate(AxeRev, P.X());
1460           if ( D.IsOpposite(gp::DX2d(),Precision::Angular())) {
1461             gp_Ax2 Ax = myCirc.Position();
1462             Ax.SetDirection(Ax.Direction().Reversed());
1463             myCirc.SetPosition(Ax);
1464           }
1465         }
1466         else if ( STy == GeomAbs_Cone) {
1467           myType = GeomAbs_Circle;
1468           gp_Cone  Cone = mySurface->Cone();
1469           gp_Pnt2d P    = myCurve->Line().Location();
1470           gp_Ax3   Axis = Cone.Position();
1471           myCirc        = ElSLib::ConeVIso(Axis,
1472                                            Cone.RefRadius(),
1473                                            Cone.SemiAngle(),
1474                                            P.Y());
1475           gp_Dir DRev = Axis.XDirection().Crossed(Axis.YDirection());
1476           gp_Ax1 AxeRev(Axis.Location(), DRev);
1477           myCirc.Rotate(AxeRev, P.X());
1478           if ( D.IsOpposite(gp::DX2d(),Precision::Angular())) {
1479             gp_Ax2 Ax = myCirc.Position();
1480             Ax.SetDirection(Ax.Direction().Reversed());
1481             myCirc.SetPosition(Ax);
1482           }
1483         }
1484         else if ( STy == GeomAbs_Torus) {
1485           myType = GeomAbs_Circle;
1486           gp_Torus Tore = mySurface->Torus();
1487           gp_Pnt2d P    = myCurve->Line().Location();
1488           gp_Ax3   Axis = Tore.Position();
1489           myCirc        = ElSLib::TorusVIso(Axis,
1490                                             Tore.MajorRadius(),
1491                                             Tore.MinorRadius(),
1492                                             P.Y());
1493           gp_Dir DRev = Axis.XDirection().Crossed(Axis.YDirection());
1494           gp_Ax1 AxeRev(Axis.Location(), DRev);
1495           myCirc.Rotate(AxeRev, P.X());
1496           if ( D.IsOpposite(gp::DX2d(),Precision::Angular())) {
1497             gp_Ax2 Ax = myCirc.Position();
1498             Ax.SetDirection(Ax.Direction().Reversed());
1499             myCirc.SetPosition(Ax);
1500           }
1501         }
1502       }
1503       else if ( D.IsParallel(gp::DY2d(),Precision::Angular())) { // Iso U.
1504         if ( STy == GeomAbs_Sphere) {
1505           myType = GeomAbs_Circle;
1506           gp_Sphere Sph  = mySurface->Sphere();
1507           gp_Pnt2d  P    = myCurve->Line().Location();
1508           gp_Ax3    Axis = Sph.Position();
1509           // calcul de l'iso 0.
1510           myCirc         = ElSLib::SphereUIso(Axis, Sph.Radius(),0.);
1511
1512           // mise a sameparameter (rotation du cercle - decalage du Y)
1513           gp_Dir DRev = Axis.XDirection().Crossed(Axis. Direction());
1514           gp_Ax1 AxeRev(Axis.Location(),DRev);
1515           myCirc.Rotate(AxeRev, P.Y());
1516
1517           // transformation en iso U ( = P.X())
1518           DRev = Axis.XDirection().Crossed(Axis.YDirection());
1519           AxeRev = gp_Ax1(Axis.Location(), DRev);
1520           myCirc.Rotate(AxeRev, P.X());
1521           
1522           if ( D.IsOpposite(gp::DY2d(),Precision::Angular())) {
1523             gp_Ax2 Ax = myCirc.Position();
1524             Ax.SetDirection(Ax.Direction().Reversed());
1525             myCirc.SetPosition(Ax);
1526           }
1527         }
1528         else if ( STy == GeomAbs_Cylinder) {
1529           myType = GeomAbs_Line;
1530           gp_Cylinder Cyl = mySurface->Cylinder();
1531           gp_Pnt2d    P   = myCurve->Line().Location();
1532           myLin           = ElSLib::CylinderUIso(Cyl.Position(),
1533                                                  Cyl.Radius(),
1534                                                  P.X());
1535           gp_Vec Tr(myLin.Direction());
1536           Tr.Multiply(P.Y());
1537           myLin.Translate(Tr);
1538           if ( D.IsOpposite(gp::DY2d(),Precision::Angular()))
1539             myLin.Reverse();
1540         }
1541         else if ( STy == GeomAbs_Cone) {
1542           myType = GeomAbs_Line;
1543           gp_Cone  Cone = mySurface->Cone();
1544           gp_Pnt2d P    = myCurve->Line().Location();
1545           myLin         = ElSLib::ConeUIso(Cone.Position(),
1546                                            Cone.RefRadius(),
1547                                            Cone.SemiAngle(),
1548                                            P.X());
1549           gp_Vec Tr(myLin.Direction());
1550           Tr.Multiply(P.Y());
1551           myLin.Translate(Tr);    
1552           if ( D.IsOpposite(gp::DY2d(),Precision::Angular()))
1553             myLin.Reverse();
1554         }
1555         else if ( STy == GeomAbs_Torus) {
1556           myType = GeomAbs_Circle;
1557           gp_Torus Tore = mySurface->Torus();
1558           gp_Pnt2d P    = myCurve->Line().Location();
1559           gp_Ax3   Axis = Tore.Position();
1560           myCirc        = ElSLib::TorusUIso(Axis,
1561                                             Tore.MajorRadius(),
1562                                             Tore.MinorRadius(),
1563                                             P.X());
1564           myCirc.Rotate(myCirc.Axis(),P.Y());
1565           
1566           if ( D.IsOpposite(gp::DY2d(),Precision::Angular())) {
1567             gp_Ax2 Ax = myCirc.Position();
1568             Ax.SetDirection(Ax.Direction().Reversed());
1569             myCirc.SetPosition(Ax);
1570           }
1571         }
1572       }
1573     }
1574   }
1575 }
1576 //=======================================================================
1577 //function :EvalFirstLastSurf 
1578 //purpose  : 
1579 //=======================================================================
1580
1581 void Adaptor3d_CurveOnSurface::EvalFirstLastSurf()
1582
1583   Standard_Real FirstPar,LastPar;
1584   gp_Pnt2d UV, LeftBot, RightTop;
1585   gp_Vec2d DUV;
1586   Standard_Real Tol= Precision::PConfusion()/10;
1587   Standard_Boolean Ok = Standard_True;
1588   
1589   
1590   FirstPar=myCurve->FirstParameter();
1591   myCurve->D1(FirstPar,UV,DUV);
1592
1593   if(DUV.Magnitude() <= Tol) Ok = Standard_False;
1594
1595   if(Ok) {
1596
1597     switch(mySurface->GetType()) {
1598     case GeomAbs_BSplineSurface :
1599       LocatePart(UV,DUV,mySurface,LeftBot,RightTop);
1600       break;
1601     case GeomAbs_SurfaceOfRevolution :
1602     case  GeomAbs_SurfaceOfExtrusion :  
1603       Ok = LocatePart_RevExt(UV,DUV,mySurface,LeftBot,RightTop);
1604       break;
1605     case GeomAbs_OffsetSurface :
1606       Ok = LocatePart_Offset(UV,DUV,mySurface,LeftBot,RightTop);
1607       break;
1608       default :
1609         Standard_NotImplemented::Raise("Adaptor3d_CurveOnSurface::EvalFirstLastSurf");
1610       break;
1611     }
1612   }
1613
1614   if (Ok) {
1615
1616     CompareBounds(LeftBot,RightTop); //SVV
1617     
1618     myFirstSurf = mySurface->UTrim(LeftBot.X(),RightTop.X(),Tol);
1619     myFirstSurf = myFirstSurf->VTrim(LeftBot.Y(),RightTop.Y(),Tol);
1620
1621   }
1622   else {
1623     myFirstSurf = mySurface;
1624   }
1625   
1626   LastPar=myCurve->LastParameter();
1627   Ok = Standard_True;
1628   myCurve->D1(LastPar,UV,DUV);
1629   DUV.Reverse(); //We want the other part
1630
1631   if(DUV.Magnitude() <= Tol) Ok = Standard_False;
1632
1633   if(Ok) {
1634
1635     switch(mySurface->GetType()) {
1636     case GeomAbs_BSplineSurface :
1637       LocatePart(UV,DUV,mySurface,LeftBot,RightTop);
1638       break;
1639     case GeomAbs_SurfaceOfRevolution :
1640     case  GeomAbs_SurfaceOfExtrusion :  
1641       Ok = LocatePart_RevExt(UV,DUV,mySurface,LeftBot,RightTop);
1642       break;
1643     case GeomAbs_OffsetSurface :
1644       Ok = LocatePart_Offset(UV,DUV,mySurface,LeftBot,RightTop);
1645       break;
1646       default :
1647         Standard_NotImplemented::Raise("Adaptor3d_CurveOnSurface::EvalFirstLastSurf");
1648       break;
1649     }
1650   }
1651
1652   if (Ok) {
1653
1654     CompareBounds(LeftBot, RightTop); //SVV
1655
1656     myLastSurf = mySurface->UTrim(LeftBot.X(),RightTop.X(),Tol);
1657     myLastSurf = myLastSurf->VTrim(LeftBot.Y(),RightTop.Y(),Tol);
1658
1659   }
1660   else {
1661     myLastSurf = mySurface;
1662   }
1663 }
1664
1665 //=======================================================================
1666 //function :LocatePart_RevExt
1667 //purpose  : processes Knots 
1668 //=======================================================================
1669
1670 Standard_Boolean Adaptor3d_CurveOnSurface::LocatePart_RevExt(const gp_Pnt2d& UV, 
1671                                                            const gp_Vec2d& DUV,
1672                                                            const Handle(Adaptor3d_HSurface)& S,
1673                                                            gp_Pnt2d& LeftBot, 
1674                                                            gp_Pnt2d& RightTop) const
1675
1676   Handle(Adaptor3d_HCurve) AHC = S->BasisCurve();
1677
1678   if (AHC->GetType() == GeomAbs_BSplineCurve) {
1679     Handle( Geom_BSplineCurve) BSplC;
1680     BSplC = AHC->BSpline();
1681  
1682     if((S->GetType())==GeomAbs_SurfaceOfExtrusion) { 
1683       Locate1Coord(1,UV,DUV,BSplC,LeftBot,RightTop);
1684       Locate2Coord(2,UV,DUV,S->FirstVParameter(),S->LastVParameter(),LeftBot,RightTop); 
1685     }
1686     else if((S->GetType())==GeomAbs_SurfaceOfRevolution)  { 
1687       Locate1Coord(2,UV,DUV,BSplC,LeftBot,RightTop);
1688       Locate2Coord(1,UV,DUV,S->FirstUParameter(),S->LastUParameter(),LeftBot,RightTop);
1689     } 
1690
1691     Standard_Real u1,u2,v1,v2;
1692     ReverseParam(LeftBot.X(),RightTop.X(),u1,u2);
1693     LeftBot.SetX(u1);
1694     RightTop.SetX(u2);
1695     ReverseParam(LeftBot.Y(),RightTop.Y(),v1,v2);
1696     LeftBot.SetY(v1);
1697     RightTop.SetY(v2);
1698     return Standard_True;
1699   }
1700   return Standard_False;
1701 }
1702
1703 //=======================================================================
1704 //function :LocatePart_OffsetSurface
1705 //purpose  : 
1706 //=======================================================================
1707
1708 Standard_Boolean Adaptor3d_CurveOnSurface::
1709    LocatePart_Offset(const gp_Pnt2d& UV, const gp_Vec2d& DUV,
1710                      const Handle(Adaptor3d_HSurface)& S,
1711                      gp_Pnt2d& LeftBot, gp_Pnt2d& RightTop) const
1712 {
1713   Standard_Boolean Ok = Standard_True;
1714   Handle( Adaptor3d_HSurface) AHS;
1715   Handle( Geom_BSplineSurface) BSplS;
1716   AHS = S->BasisSurface();
1717   GeomAbs_SurfaceType  BasisSType = AHS->GetType();
1718   switch(BasisSType) {
1719   case GeomAbs_SurfaceOfRevolution:
1720   case GeomAbs_SurfaceOfExtrusion :
1721     Ok = LocatePart_RevExt(UV,DUV,AHS,LeftBot,RightTop);
1722     break;
1723       
1724   case GeomAbs_BSplineSurface:
1725     LocatePart(UV,DUV,AHS,LeftBot,RightTop);
1726     break;
1727
1728     default : 
1729       Ok=Standard_False;
1730   }
1731   return Ok;
1732 }
1733
1734 //=======================================================================
1735 //function :LocatePart
1736 //purpose  : for BSplineSurface
1737 //=======================================================================
1738
1739 void Adaptor3d_CurveOnSurface::LocatePart(const gp_Pnt2d& UV, const gp_Vec2d& DUV,
1740                                         const Handle(Adaptor3d_HSurface)& S,
1741                                         gp_Pnt2d& LeftBot, gp_Pnt2d& RightTop) const
1742 {
1743   Handle( Geom_BSplineSurface) BSplS;
1744   BSplS = S->BSpline();
1745   Standard_Boolean  DUIsNull=Standard_False,
1746   DVIsNull=Standard_False;
1747       
1748   Locate1Coord(1,UV,DUV,BSplS,DUIsNull,LeftBot,RightTop);
1749   Locate1Coord(2,UV,DUV,BSplS,DVIsNull,LeftBot,RightTop); 
1750   
1751   if((DUIsNull==Standard_True)&&(DVIsNull==Standard_False)) {
1752     TColStd_Array1OfReal ArrU(1,BSplS->NbUKnots());
1753     BSplS->UKnots(ArrU);
1754     Locate2Coord(1,UV,DUV,BSplS,ArrU,LeftBot,RightTop); 
1755   }
1756   else if((DVIsNull==Standard_True)&&(DUIsNull==Standard_False)) {
1757     TColStd_Array1OfReal ArrV(1,BSplS->NbVKnots());
1758     BSplS->VKnots(ArrV);
1759     Locate2Coord(2,UV,DUV,BSplS,ArrV,LeftBot,RightTop);
1760     }
1761 }
1762
1763
1764