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