src/Adaptor3d/Adaptor3d_CurveOnSurface.cxx

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