src/Adaptor3d/Adaptor3d_CurveOnSurface.cxx

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