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