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