src/Adaptor3d/Adaptor3d_CurveOnSurface.cxx

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