src/Adaptor3d/Adaptor3d_CurveOnSurface.cxx

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