src/Adaptor3d/Adaptor3d_CurveOnSurface.cxx

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