src/Geom/Geom_OffsetSurface.cxx

   1 // Created on: 1991-06-25
   2 // Created by: JCV
   3 // Copyright (c) 1991-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 // Modified     04/10/96 : JCT : derivee des surfaces offset utilisation de
  18 //                               CSLib
  19 // Modified     15/11/96 : JPI : ajout equivalent surface pour les surfaces canoniques et modif des methodes D0 D1, ... UIso,VIso
  20 // Modified     18/11/96 : JPI : inversion de l'offsetValue dans UReverse et Vreverse
  21
  22 #include <AdvApprox_ApproxAFunction.hxx>
  23 #include <BSplCLib.hxx>
  24 #include <BSplSLib.hxx>
  25 #include <Convert_GridPolynomialToPoles.hxx>
  26 #include <CSLib.hxx>
  27 #include <Geom_BezierSurface.hxx>
  28 #include <Geom_BSplineCurve.hxx>
  29 #include <Geom_BSplineSurface.hxx>
  30 #include <Geom_Circle.hxx>
  31 #include <Geom_ConicalSurface.hxx>
  32 #include <Geom_Curve.hxx>
  33 #include <Geom_CylindricalSurface.hxx>
  34 #include <Geom_ElementarySurface.hxx>
  35 #include <Geom_Ellipse.hxx>
  36 #include <Geom_Geometry.hxx>
  37 #include <Geom_OffsetCurve.hxx>
  38 #include <Geom_OffsetSurface.hxx>
  39 #include <Geom_Plane.hxx>
  40 #include <Geom_RectangularTrimmedSurface.hxx>
  41 #include <Geom_SphericalSurface.hxx>
  42 #include <Geom_Surface.hxx>
  43 #include <Geom_SurfaceOfLinearExtrusion.hxx>
  44 #include <Geom_SurfaceOfRevolution.hxx>
  45 #include <Geom_ToroidalSurface.hxx>
  46 #include <Geom_TrimmedCurve.hxx>
  47 #include <Geom_UndefinedDerivative.hxx>
  48 #include <Geom_UndefinedValue.hxx>
  49 #include <GeomAbs_CurveType.hxx>
  50 #include <GeomAbs_IsoType.hxx>
  51 #include <GeomAbs_Shape.hxx>
  52 #include <gp.hxx>
  53 #include <gp_Dir.hxx>
  54 #include <gp_GTrsf2d.hxx>
  55 #include <gp_Pnt.hxx>
  56 #include <gp_Trsf.hxx>
  57 #include <gp_Vec.hxx>
  58 #include <gp_XYZ.hxx>
  59 #include <Precision.hxx>
  60 #include <Standard_ConstructionError.hxx>
  61 #include <Standard_NoSuchObject.hxx>
  62 #include <Standard_NotImplemented.hxx>
  63 #include <Standard_RangeError.hxx>
  64 #include <Standard_Type.hxx>
  65 #include <TColgp_Array1OfPnt.hxx>
  66 #include <TColgp_Array2OfVec.hxx>
  67 #include <TColgp_HArray2OfPnt.hxx>
  68 #include <TColStd_Array1OfInteger.hxx>
  69 #include <TColStd_Array1OfReal.hxx>
  70 #include <TColStd_HArray1OfInteger.hxx>
  71 #include <TColStd_HArray1OfReal.hxx>
  72 #include <TColStd_HArray2OfInteger.hxx>
  73
  74 static const Standard_Real MyAngularToleranceForG1 = Precision::Angular();
  75
  76 //=======================================================================
  77 //function : derivatives
  78 //purpose  :
  79 //=======================================================================
  80
  81 static void derivatives(Standard_Integer MaxOrder,
  82   Standard_Integer MinOrder,
  83   const Standard_Real U,
  84   const Standard_Real V,
  85   const Handle(Geom_Surface)& basisSurf,
  86   const Standard_Integer Nu,
  87   const Standard_Integer Nv,
  88   const Standard_Boolean AlongU,
  89   const Standard_Boolean AlongV,
  90   const Handle(Geom_BSplineSurface)& L,
  91   TColgp_Array2OfVec& DerNUV,
  92   TColgp_Array2OfVec& DerSurf)
  93 {
  94   Standard_Integer i,j;
  95   gp_Pnt P;
  96   gp_Vec DL1U, DL1V, DL2U , DL2V , DL2UV ,DL3U, DL3UUV, DL3UVV, DL3V;
  97
  98   if (AlongU || AlongV)
  99   {
 100     MaxOrder=0;
 101     TColgp_Array2OfVec DerSurfL(0,MaxOrder+Nu+1,0,MaxOrder+Nv+1);
 102     switch (MinOrder)
 103     {
 104     case 1 :
 105       L->D1(U, V, P, DL1U, DL1V);
 106       DerSurfL.SetValue(1, 0, DL1U);
 107       DerSurfL.SetValue(0, 1, DL1V);
 108       break;
 109     case 2 :
 110       L->D2(U, V, P, DL1U, DL1V, DL2U , DL2V , DL2UV);
 111       DerSurfL.SetValue(1, 0, DL1U);
 112       DerSurfL.SetValue(0, 1, DL1V);
 113       DerSurfL.SetValue(1, 1, DL2UV);
 114       DerSurfL.SetValue(2, 0, DL2U);
 115       DerSurfL.SetValue(0, 2, DL2V);
 116       break;
 117     case 3 :
 118       L->D3(U, V, P, DL1U, DL1V, DL2U , DL2V , DL2UV ,DL3U, DL3V ,DL3UUV, DL3UVV);
 119       DerSurfL.SetValue(1, 0, DL1U);
 120       DerSurfL.SetValue(0, 1, DL1V);
 121       DerSurfL.SetValue(1, 1, DL2UV);
 122       DerSurfL.SetValue(2, 0, DL2U);
 123       DerSurfL.SetValue(0, 2, DL2V);
 124       DerSurfL.SetValue(3, 0, DL3U);
 125       DerSurfL.SetValue(2, 1, DL3UUV);
 126       DerSurfL.SetValue(1, 2, DL3UVV);
 127       DerSurfL.SetValue(0, 3, DL3V);
 128       break;
 129     default:
 130       break;
 131     }
 132
 133     if(Nu <= Nv)
 134     {
 135       for(i=0;i<=MaxOrder+1+Nu;i++)
 136         for(j=i;j<=MaxOrder+Nv+1;j++)
 137           if(i+j> MinOrder)
 138           {
 139             DerSurfL.SetValue(i,j,L->DN(U,V,i,j));
 140             DerSurf.SetValue(i,j,basisSurf->DN(U,V,i,j));
 141             if (i!=j && j <= Nu+1)
 142             {
 143               DerSurf.SetValue(j,i,basisSurf->DN(U,V,j,i));
 144               DerSurfL.SetValue(j,i,L->DN(U,V,j,i));
 145             }
 146           }
 147     }
 148     else
 149     {
 150       for(j=0;j<=MaxOrder+1+Nv;j++)
 151         for(i=j;i<=MaxOrder+Nu+1;i++)
 152           if(i+j> MinOrder)
 153           {
 154             DerSurfL.SetValue(i,j,L->DN(U,V,i,j));
 155             DerSurf.SetValue(i,j,basisSurf->DN(U,V,i,j));
 156             if (i!=j && i <= Nv+1)
 157             {
 158               DerSurf.SetValue(j,i,basisSurf->DN(U,V,j,i));
 159               DerSurfL.SetValue(j,i,L->DN(U,V,j,i));
 160             }
 161           }
 162     }
 163     for(i=0;i<=MaxOrder+Nu;i++)
 164       for(j=0;j<=MaxOrder+Nv;j++)
 165       {
 166         if (AlongU)
 167           DerNUV.SetValue(i,j,CSLib::DNNUV(i,j,DerSurfL,DerSurf));
 168         if (AlongV)
 169           DerNUV.SetValue(i,j,CSLib::DNNUV(i,j,DerSurf,DerSurfL));
 170       }
 171
 172   }
 173   else
 174   {
 175
 176     for(i=0;i<=MaxOrder+Nu+1;i++)
 177       for(j=i;j<=MaxOrder+Nv+1;j++)
 178         if(i+j>MinOrder)
 179         {
 180
 181           DerSurf.SetValue(i,j,basisSurf->DN(U,V,i,j));
 182           if (i!=j) DerSurf.SetValue(j,i,basisSurf->DN(U,V,j,i));
 183         }
 184         for(i=0;i<=MaxOrder+Nu;i++)
 185           for(j=0;j<=MaxOrder+Nv;j++)
 186           {
 187             DerNUV.SetValue(i,j,CSLib::DNNUV(i,j,DerSurf));
 188           }
 189   }
 190
 191 }
 192
 193 //=======================================================================
 194 //function : Copy
 195 //purpose  :
 196 //=======================================================================
 197
 198 Handle(Geom_Geometry) Geom_OffsetSurface::Copy () const
 199 {
 200   Handle(Geom_OffsetSurface) S(new Geom_OffsetSurface(basisSurf, offsetValue, Standard_True));
 201   return S;
 202 }
 203
 204 //=======================================================================
 205 //function : Geom_OffsetSurface
 206 //purpose  : Basis surface cannot be an Offset surface or trimmed from
 207 //            offset surface.
 208 //=======================================================================
 209
 210 Geom_OffsetSurface::Geom_OffsetSurface (const Handle(Geom_Surface)& theSurf,
 211   const Standard_Real theOffset,
 212   const Standard_Boolean isNotCheckC0)
 213   : offsetValue (theOffset)
 214 {
 215   SetBasisSurface(theSurf, isNotCheckC0);
 216 }
 217
 218 //=======================================================================
 219 //function : SetBasisSurface
 220 //purpose  :
 221 //=======================================================================
 222
 223 void Geom_OffsetSurface::SetBasisSurface (const Handle(Geom_Surface)& S,
 224   const Standard_Boolean isNotCheckC0)
 225 {
 226   Standard_Real aUf, aUl, aVf, aVl;
 227   S->Bounds(aUf, aUl, aVf, aVl);
 228
 229   Handle(Geom_Surface) aCheckingSurf = Handle(Geom_Surface)::DownCast(S->Copy());
 230   Standard_Boolean isTrimmed = Standard_False;
 231
 232   while(aCheckingSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
 233         aCheckingSurf->IsKind(STANDARD_TYPE(Geom_OffsetSurface)))
 234   {
 235     if (aCheckingSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
 236     {
 237       Handle(Geom_RectangularTrimmedSurface) aTrimS =
 238         Handle(Geom_RectangularTrimmedSurface)::DownCast(aCheckingSurf);
 239       aCheckingSurf = aTrimS->BasisSurface();
 240       isTrimmed = Standard_True;
 241     }
 242
 243     if (aCheckingSurf->IsKind(STANDARD_TYPE(Geom_OffsetSurface)))
 244     {
 245       Handle(Geom_OffsetSurface) aOS =
 246         Handle(Geom_OffsetSurface)::DownCast(aCheckingSurf);
 247       aCheckingSurf = aOS->BasisSurface();
 248       offsetValue += aOS->Offset();
 249     }
 250   }
 251
 252   myBasisSurfContinuity = aCheckingSurf->Continuity();
 253
 254   Standard_Boolean isC0 = !isNotCheckC0 && (myBasisSurfContinuity == GeomAbs_C0);
 255
 256   // Basis surface must be at least C1
 257   if (isC0)
 258   {
 259     Handle(Geom_Curve) aCurve;
 260
 261     if (aCheckingSurf->IsKind(STANDARD_TYPE(Geom_SurfaceOfRevolution)))
 262     {
 263       Handle(Geom_SurfaceOfRevolution) aRevSurf = Handle(Geom_SurfaceOfRevolution)::DownCast(aCheckingSurf);
 264       aCurve = aRevSurf->BasisCurve();
 265     }
 266     else if (aCheckingSurf->IsKind(STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion)))
 267     {
 268       Handle(Geom_SurfaceOfLinearExtrusion) aLESurf = Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(aCheckingSurf);
 269       aCurve = aLESurf->BasisCurve();
 270     }
 271
 272     if(!aCurve.IsNull())
 273     {
 274       while(aCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)) ||
 275             aCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve)))
 276       {
 277         if (aCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)))
 278         {
 279           Handle(Geom_TrimmedCurve) aTrimC =
 280             Handle(Geom_TrimmedCurve)::DownCast(aCurve);
 281           aCurve = aTrimC->BasisCurve();
 282         }
 283
 284         if (aCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve)))
 285         {
 286           Handle(Geom_OffsetCurve) aOC =
 287             Handle(Geom_OffsetCurve)::DownCast(aCurve);
 288           aCurve = aOC->BasisCurve();
 289         }
 290       }
 291     }
 292
 293     const Standard_Real aUIsoPar = (aUf + aUl)/2.0, aVIsoPar = (aVf + aVl)/2.0;
 294     Standard_Boolean isUG1 = Standard_False, isVG1 = Standard_False;
 295
 296     const Handle(Geom_Curve) aCurv1 = aCurve.IsNull() ? aCheckingSurf->UIso(aUIsoPar) : aCurve;
 297     const Handle(Geom_Curve) aCurv2 = aCheckingSurf->VIso(aVIsoPar);
 298     isUG1 = !aCurv1->IsKind(STANDARD_TYPE(Geom_BSplineCurve));
 299     isVG1 = !aCurv2->IsKind(STANDARD_TYPE(Geom_BSplineCurve));
 300
 301     if(!isUG1)
 302     {
 303       Handle(Geom_BSplineCurve) aBC = Handle(Geom_BSplineCurve)::DownCast(aCurv1);
 304       isUG1 = aBC->IsG1(aVf, aVl, MyAngularToleranceForG1);
 305     }
 306     //
 307     if(!isVG1)
 308     {
 309       Handle(Geom_BSplineCurve) aBC = Handle(Geom_BSplineCurve)::DownCast(aCurv2);
 310       isVG1 = aBC->IsG1(aUf, aUl, MyAngularToleranceForG1);
 311     }
 312     //
 313     if(isUG1 && isVG1)
 314     {
 315       myBasisSurfContinuity = GeomAbs_G1;
 316       isC0 = Standard_False;
 317     }
 318
 319     // Raise exception if still C0
 320     if (isC0)
 321       Standard_ConstructionError::Raise("Offset with no C1 Surface");
 322   }
 323
 324   if(isTrimmed)
 325   {
 326     basisSurf =
 327       new Geom_RectangularTrimmedSurface(aCheckingSurf, aUf, aUl, aVf, aVl);
 328   }
 329   else
 330   {
 331     basisSurf = aCheckingSurf;
 332   }
 333
 334   equivSurf = Surface();
 335   //
 336   // Tolerance en dur pour l'instant ,mais on devrait la proposer dans le constructeur
 337   // et la mettre en champ, on pourrait utiliser par exemple pour l'extraction d'iso
 338   // et aussi pour les singularite. Pour les surfaces osculatrices, on l'utilise pour
 339   // detecter si une iso est degeneree.
 340   const Standard_Real Tol = Precision::Confusion(); //0.0001;
 341   myOscSurf.Init(basisSurf,Tol);
 342 }
 343
 344 //=======================================================================
 345 //function : SetOffsetValue
 346 //purpose  :
 347 //=======================================================================
 348
 349 void Geom_OffsetSurface::SetOffsetValue (const Standard_Real D)
 350 {
 351   offsetValue = D;
 352   equivSurf = Surface();
 353 }
 354
 355 //=======================================================================
 356 //function : UReverse
 357 //purpose  :
 358 //=======================================================================
 359
 360 void Geom_OffsetSurface::UReverse ()
 361 {
 362   basisSurf->UReverse();
 363   offsetValue = -offsetValue;
 364   if(!equivSurf.IsNull()) equivSurf->UReverse();
 365 }
 366
 367 //=======================================================================
 368 //function : UReversedParameter
 369 //purpose  :
 370 //=======================================================================
 371
 372 Standard_Real Geom_OffsetSurface::UReversedParameter(const Standard_Real U) const
 373 {
 374   return basisSurf->UReversedParameter(U);
 375 }
 376
 377 //=======================================================================
 378 //function : VReverse
 379 //purpose  :
 380 //=======================================================================
 381
 382 void Geom_OffsetSurface::VReverse ()
 383 {
 384   basisSurf->VReverse();
 385   offsetValue = -offsetValue;
 386   if(!equivSurf.IsNull()) equivSurf->VReverse();
 387 }
 388
 389 //=======================================================================
 390 //function : VReversedParameter
 391 //purpose  :
 392 //=======================================================================
 393
 394 Standard_Real Geom_OffsetSurface::VReversedParameter(const Standard_Real V) const
 395 {
 396   return basisSurf->VReversedParameter(V);
 397 }
 398
 399 //=======================================================================
 400 //function : Bounds
 401 //purpose  :
 402 //=======================================================================
 403
 404 void Geom_OffsetSurface::Bounds (Standard_Real& U1, Standard_Real& U2,
 405                                  Standard_Real& V1, Standard_Real& V2) const
 406 {
 407   basisSurf->Bounds (U1, U2 ,V1, V2);
 408 }
 409
 410 //=======================================================================
 411 //function : Continuity
 412 //purpose  :
 413 //=======================================================================
 414
 415 GeomAbs_Shape Geom_OffsetSurface::Continuity () const
 416 {
 417   switch (myBasisSurfContinuity) {
 418     case GeomAbs_C2 : return GeomAbs_C1;
 419     case GeomAbs_C3 : return GeomAbs_C2;
 420     case GeomAbs_CN : return GeomAbs_CN;
 421     default : break;
 422   }
 423   return GeomAbs_C0;
 424 }
 425
 426 //=======================================================================
 427 //function : D0
 428 //purpose  :
 429 //=======================================================================
 430
 431 void Geom_OffsetSurface::D0 (const Standard_Real U, const Standard_Real V, gp_Pnt& P) const
 432 {
 433   gp_Vec D1U, D1V;
 434 #ifdef CHECK
 435   if (myBasisSurfContinuity == GeomAbs_C0)
 436     Geom_UndefinedValue::Raise();
 437 #endif
 438   if (equivSurf.IsNull()){
 439     basisSurf->D1(U, V, P, D1U, D1V);
 440     SetD0(U,V,P,D1U,D1V);
 441   }
 442   else
 443     equivSurf->D0(U,V,P);
 444 }
 445
 446 //=======================================================================
 447 //function : D1
 448 //purpose  :
 449 //=======================================================================
 450
 451 void Geom_OffsetSurface::D1 (const Standard_Real U, const Standard_Real V,
 452   gp_Pnt& P,
 453   gp_Vec& D1U, gp_Vec& D1V) const
 454 {
 455 #ifdef CHECK
 456   if (myBasisSurfContinuity == GeomAbs_C0 ||
 457       myBasisSurfContinuity == GeomAbs_C1)
 458     Geom_UndefinedDerivative::Raise();
 459 #endif
 460   if (equivSurf.IsNull())
 461   {
 462     gp_Vec d2u, d2v, d2uv;
 463     basisSurf->D2(U, V, P, D1U, D1V, d2u, d2v, d2uv);
 464     SetD1(U,V,P,D1U,D1V,d2u,d2v,d2uv);
 465   }
 466   else
 467     equivSurf->D1(U,V,P,D1U,D1V);
 468 }
 469
 470 //=======================================================================
 471 //function : D2
 472 //purpose  :
 473 //=======================================================================
 474
 475 void Geom_OffsetSurface::D2 (const Standard_Real U, const Standard_Real V,
 476   gp_Pnt& P,
 477   gp_Vec& D1U, gp_Vec& D1V,
 478   gp_Vec& D2U, gp_Vec& D2V, gp_Vec& D2UV) const
 479 {
 480 #ifdef CHECK
 481   if (myBasisSurfContinuity == GeomAbs_C0 ||
 482       myBasisSurfContinuity == GeomAbs_C1 ||
 483       myBasisSurfContinuity == GeomAbs_C2)
 484     Geom_UndefinedDerivative::Raise();
 485 #endif
 486   if (equivSurf.IsNull())
 487   {
 488     gp_Vec d3u, d3uuv, d3uvv, d3v;
 489     basisSurf->D3(U, V, P, D1U, D1V, D2U, D2V, D2UV, d3u, d3v, d3uuv, d3uvv);
 490     SetD2(U,V,P,D1U,D1V,D2U,D2V,D2UV,d3u,d3v,d3uuv,d3uvv);
 491   }
 492   else
 493     equivSurf->D2(U,V,P,D1U,D1V,D2U,D2V,D2UV);
 494 }
 495
 496 //=======================================================================
 497 //function : D3
 498 //purpose  :
 499 //=======================================================================
 500
 501 void Geom_OffsetSurface::D3 (const Standard_Real U, const Standard_Real V,
 502   gp_Pnt& P,
 503   gp_Vec& D1U, gp_Vec& D1V,
 504   gp_Vec& D2U, gp_Vec& D2V, gp_Vec& D2UV,
 505   gp_Vec& D3U, gp_Vec& D3V, gp_Vec& D3UUV, gp_Vec& D3UVV) const
 506 {
 507 #ifdef CHECK
 508   if (!(basisSurf->IsCNu (4) && basisSurf->IsCNv (4))) {
 509     Geom_UndefinedDerivative::Raise();
 510   }
 511 #endif
 512   if (equivSurf.IsNull())
 513   {
 514     basisSurf->D3(U, V, P, D1U, D1V, D2U, D2V, D2UV, D3U, D3V, D3UUV, D3UVV);
 515     SetD3(U, V, P, D1U, D1V, D2U, D2V, D2UV, D3U, D3V, D3UUV, D3UVV);
 516   }
 517   else
 518     equivSurf->D3(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 519 }
 520
 521 //=======================================================================
 522 //function : DN
 523 //purpose  :
 524 //=======================================================================
 525
 526 gp_Vec Geom_OffsetSurface::DN (const Standard_Real    U, const Standard_Real    V,
 527   const Standard_Integer Nu, const Standard_Integer Nv) const
 528 {
 529   Standard_RangeError_Raise_if (Nu < 0 || Nv < 0 || Nu + Nv < 1, " ");
 530 #ifdef CHECK
 531   if (!(basisSurf->IsCNu (Nu) && basisSurf->IsCNv (Nv))) {
 532     Geom_UndefinedDerivative::Raise();
 533   }
 534 #endif
 535   gp_Vec D(0,0,0);
 536
 537   if (equivSurf.IsNull())
 538   {
 539     gp_Pnt P;
 540     gp_Vec D1U,D1V;
 541     basisSurf->D1 (U, V, P, D1U, D1V);
 542
 543     D = SetDN(U,V,Nu,Nv,D1U,D1V);
 544   }
 545   else
 546     D = equivSurf->DN(U,V,Nu,Nv);
 547   return D;
 548 }
 549
 550 //=======================================================================
 551 //function : D1
 552 //purpose  :
 553 //=======================================================================
 554
 555 void Geom_OffsetSurface::D1
 556   (const Standard_Real U,  const Standard_Real V,
 557    gp_Pnt& P,        gp_Pnt& Pbasis,
 558    gp_Vec& D1U,      gp_Vec& D1V,      gp_Vec& D1Ubasis, gp_Vec& D1Vbasis,
 559    gp_Vec& D2Ubasis, gp_Vec& D2Vbasis, gp_Vec& D2UVbasis) const
 560 {
 561   const GeomAbs_Shape basisCont = basisSurf->Continuity();
 562   if (basisCont == GeomAbs_C0 || basisCont == GeomAbs_C1)
 563     Geom_UndefinedDerivative::Raise();
 564   basisSurf->D2 (U, V, Pbasis, D1Ubasis, D1Vbasis, D2Ubasis, D2Vbasis, D2UVbasis);
 565   gp_Vec Ndir = D1Ubasis.Crossed (D1Vbasis);
 566   const Standard_Real R2  = Ndir.SquareMagnitude();
 567   const Standard_Real R   = Sqrt (R2);
 568   const Standard_Real R3  = R * R2;
 569   gp_Vec DUNdir = D2Ubasis.Crossed (D1Vbasis);
 570   DUNdir.Add (D1Ubasis.Crossed (D2UVbasis));
 571   gp_Vec DVNdir = D2UVbasis.Crossed (D1Vbasis);
 572   DVNdir.Add (D1Ubasis.Crossed (D2Vbasis));
 573   const Standard_Real DRu = Ndir.Dot (DUNdir);
 574   const Standard_Real DRv = Ndir.Dot (DVNdir);
 575   if (R3 <= gp::Resolution()) {
 576     if (R2 <= gp::Resolution())
 577       Geom_UndefinedDerivative::Raise();
 578     DUNdir.Multiply(R);
 579     DUNdir.Subtract (Ndir.Multiplied (DRu/R));
 580     DUNdir.Multiply (offsetValue/R2);
 581     D1U = D1Ubasis.Added (DUNdir);
 582     DVNdir.Multiply(R);
 583     DVNdir.Subtract (Ndir.Multiplied (DRv/R));
 584     DVNdir.Multiply (offsetValue/R2);
 585     D1V = D1Vbasis.Added (DVNdir);
 586   }
 587   else {
 588     DUNdir.Multiply (offsetValue / R);
 589     DUNdir.Subtract (Ndir.Multiplied (offsetValue*DRu/R3));
 590     D1U = D1Ubasis.Added (DUNdir);
 591     DVNdir.Multiply (offsetValue / R);
 592     DVNdir.Subtract (Ndir.Multiplied (offsetValue*DRv/R3));
 593     D1V = D1Vbasis.Added (DVNdir);
 594   }
 595   Ndir.Multiply (offsetValue/R);
 596   P.SetXYZ ((Ndir.XYZ()).Added (Pbasis.XYZ()));
 597 }
 598
 599 //=======================================================================
 600 //function : D2
 601 //purpose  :
 602 //=======================================================================
 603
 604 void Geom_OffsetSurface::D2
 605   (const Standard_Real U, const Standard_Real V,
 606    gp_Pnt& P, gp_Pnt& Pbasis,
 607    gp_Vec& D1U, gp_Vec& D1V, gp_Vec& D2U, gp_Vec& D2V, gp_Vec& D2UV,
 608    gp_Vec& D1Ubasis, gp_Vec& D1Vbasis,
 609    gp_Vec& D2Ubasis, gp_Vec& D2Vbasis, gp_Vec& D2UVbasis,
 610    gp_Vec& D3Ubasis, gp_Vec& D3Vbasis, gp_Vec& D3UUVbasis, gp_Vec& D3UVVbasis) const
 611 {
 612   const GeomAbs_Shape basisCont = basisSurf->Continuity();
 613   if (basisCont == GeomAbs_C0 || basisCont == GeomAbs_C1 || basisCont == GeomAbs_C2)
 614     Geom_UndefinedDerivative::Raise();
 615   basisSurf->D3 (U, V, P, D1Ubasis, D1Vbasis, D2Ubasis, D2Vbasis, D2UVbasis, D3Ubasis, D3Vbasis, D3UUVbasis, D3UVVbasis);
 616   gp_Vec Ndir = D1Ubasis.Crossed (D1Vbasis);
 617   const Standard_Real R2  = Ndir.SquareMagnitude();
 618   const Standard_Real R   = Sqrt (R2);
 619   const Standard_Real R3  = R2 * R;
 620   const Standard_Real R4  = R2 * R2;
 621   const Standard_Real R5  = R3 * R2;
 622   gp_Vec DUNdir = D2Ubasis.Crossed (D1Vbasis);
 623   DUNdir.Add (D1Ubasis.Crossed (D2UVbasis));
 624   gp_Vec DVNdir = D2UVbasis.Crossed (D1Vbasis);
 625   DVNdir.Add (D1Ubasis.Crossed (D2Vbasis));
 626   const Standard_Real DRu = Ndir.Dot (DUNdir);
 627   const Standard_Real DRv = Ndir.Dot (DVNdir);
 628   gp_Vec D2UNdir = D3Ubasis.Crossed (D1Vbasis);
 629   D2UNdir.Add (D1Ubasis.Crossed (D3UUVbasis));
 630   D2UNdir.Add (2.0 * (D2Ubasis.Crossed (D2UVbasis)));
 631   gp_Vec D2VNdir = D3UVVbasis.Crossed (D1Vbasis);
 632   D2VNdir.Add (D1Ubasis.Crossed (D3Vbasis));
 633   D2VNdir.Add (2.0 * (D2UVbasis.Crossed (D2Vbasis)));
 634   gp_Vec D2UVNdir = D2UVbasis.Crossed (D1Vbasis);
 635   D2UVNdir.Add (D1Ubasis.Crossed (D3UVVbasis));
 636   D2UVNdir.Add (D2Ubasis.Crossed (D2Vbasis));
 637   const Standard_Real D2Ru = Ndir.Dot (D2UNdir) + DUNdir.Dot (DUNdir);
 638   const Standard_Real D2Rv = Ndir.Dot (D2VNdir) + DVNdir.Dot (DVNdir);
 639   const Standard_Real D2Ruv = DVNdir.Dot (DUNdir) + Ndir.Dot (D2UVNdir);
 640
 641   if (R5 <= gp::Resolution()) {
 642     //We try another computation but the stability is not very good
 643     //dixit ISG.
 644     if (R4 <= gp::Resolution()) Geom_UndefinedDerivative::Raise();
 645     // V2 = P" (U) :
 646     // Standard_Real R4 = R2 * R2;
 647     D2UNdir.Subtract (DUNdir.Multiplied (2.0 * DRu / R2));
 648     D2UNdir.Subtract (Ndir.Multiplied (D2Ru/R2));
 649     D2UNdir.Add (Ndir.Multiplied (3.0 * DRu * DRu / R4));
 650     D2UNdir.Multiply (offsetValue / R);
 651     D2U = D2Ubasis.Added (D2UNdir);
 652     D2VNdir.Subtract (DVNdir.Multiplied (2.0 * DRv / R2));
 653     D2VNdir.Subtract (Ndir.Multiplied (D2Rv/R2));
 654     D2VNdir.Add (Ndir.Multiplied (3.0 * DRv * DRv / R4));
 655     D2VNdir.Multiply (offsetValue / R);
 656     D2V = D2Vbasis.Added (D2VNdir);
 657
 658     D2UVNdir.Subtract (DUNdir.Multiplied (DRv / R2));
 659     D2UVNdir.Subtract (DVNdir.Multiplied (DRu / R2));
 660     D2UVNdir.Subtract (Ndir.Multiplied (D2Ruv / R2));
 661     D2UVNdir.Add (Ndir.Multiplied (3.0 * DRu * DRv / R4));
 662     D2UVNdir.Multiply (offsetValue / R);
 663     D2UV = D2UVbasis.Added (D2UVNdir);
 664
 665     DUNdir.Multiply(R);
 666     DUNdir.Subtract (Ndir.Multiplied (DRu/R));
 667     DUNdir.Multiply (offsetValue/R2);
 668     D1U = D1Ubasis.Added (DUNdir);
 669     DVNdir.Multiply(R);
 670     DVNdir.Subtract (Ndir.Multiplied (DRv/R));
 671     DVNdir.Multiply (offsetValue/R2);
 672     D1V = D1Vbasis.Added (DVNdir);
 673   }
 674   else {
 675     D2UNdir.Multiply (offsetValue/R);
 676     D2UNdir.Subtract (DUNdir.Multiplied (2.0 * offsetValue * DRu / R3));
 677     D2UNdir.Subtract (Ndir.Multiplied (offsetValue * D2Ru / R3));
 678     D2UNdir.Add (Ndir.Multiplied (offsetValue * 3.0 * DRu * DRu / R5));
 679     D2U = D2Ubasis.Added (D2UNdir);
 680
 681     D2VNdir.Multiply (offsetValue/R);
 682     D2VNdir.Subtract (DVNdir.Multiplied (2.0 * offsetValue * DRv / R3));
 683     D2VNdir.Subtract (Ndir.Multiplied (offsetValue * D2Rv / R3));
 684     D2VNdir.Add (Ndir.Multiplied (offsetValue * 3.0 * DRv * DRv / R5));
 685     D2V = D2Vbasis.Added (D2VNdir);
 686
 687     D2UVNdir.Multiply (offsetValue/R);
 688     D2UVNdir.Subtract (DUNdir.Multiplied (offsetValue * DRv / R3));
 689     D2UVNdir.Subtract (DVNdir.Multiplied (offsetValue * DRu / R3));
 690     D2UVNdir.Subtract (Ndir.Multiplied (offsetValue * D2Ruv / R3));
 691     D2UVNdir.Add (Ndir.Multiplied (3.0 * offsetValue * DRu * DRv / R5));
 692     D2UV = D2UVbasis.Added (D2UVNdir);
 693
 694     DUNdir.Multiply (offsetValue / R);
 695     DUNdir.Subtract (Ndir.Multiplied (offsetValue*DRu/R3));
 696     D1U = D1Ubasis.Added (DUNdir);
 697     DVNdir.Multiply (offsetValue / R);
 698     DVNdir.Subtract (Ndir.Multiplied (offsetValue*DRv/R3));
 699     D1V = D1Vbasis.Added (DVNdir);
 700   }
 701   Ndir.Multiply (offsetValue/R);
 702   P.SetXYZ ((Ndir.XYZ()).Added (Pbasis.XYZ()));
 703 }
 704
 705 //=======================================================================
 706 //function : LocalD0
 707 //purpose  :
 708 //=======================================================================
 709
 710 void Geom_OffsetSurface::LocalD0 (const Standard_Real U, const Standard_Real V,
 711   const Standard_Integer USide, const Standard_Integer VSide,
 712   gp_Pnt& P) const
 713 {
 714   if (equivSurf.IsNull()) {
 715     gp_Vec D1U, D1V;
 716     Handle(Geom_Surface) Basis = basisSurf;
 717
 718     // if Basis is Trimmed we take the basis's basis
 719     Handle(Geom_RectangularTrimmedSurface) RTS =
 720       Handle(Geom_RectangularTrimmedSurface)::DownCast(Basis);
 721     if (!RTS.IsNull())
 722       Basis = RTS->BasisSurface();
 723
 724     // BSpline case
 725     Handle(Geom_BSplineSurface) BSplS =
 726       Handle(Geom_BSplineSurface)::DownCast(Basis);
 727     if (!BSplS.IsNull()) {
 728       gp_Vec D2U,D2V,D2UV, D3U,D3V,D3UUV,D3UVV;
 729       LocateSides(U,V,USide,VSide,BSplS,1,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 730       SetD0(U,V,P,D1U,D1V);
 731       return;
 732     }
 733
 734     // Extrusion case
 735     Handle( Geom_SurfaceOfLinearExtrusion) SE =
 736       Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(Basis);
 737     if (!SE.IsNull()) {
 738       SE->LocalD1(U,V,USide,P,D1U,D1V);
 739       SetD0(U,V,P,D1U,D1V);
 740       return;
 741     }
 742
 743     // Revolution case
 744     Handle(Geom_SurfaceOfRevolution) SR =
 745       Handle(Geom_SurfaceOfRevolution)::DownCast(Basis);
 746     if (!SR.IsNull()) {
 747       SR->LocalD1(U,V,VSide,P,D1U,D1V);
 748       SetD0(U,V,P,D1U,D1V);
 749       return;
 750     }
 751
 752     // General cases
 753     basisSurf->D1(U, V, P, D1U, D1V);
 754     SetD0(U,V,P,D1U,D1V);
 755   }
 756   else
 757     equivSurf-> D0(U,V,P);
 758 }
 759
 760 //=======================================================================
 761 //function : LocalD1
 762 //purpose  :
 763 //=======================================================================
 764
 765 void Geom_OffsetSurface::LocalD1 (const Standard_Real U, const Standard_Real V,
 766   const Standard_Integer USide, const Standard_Integer VSide,
 767   gp_Pnt& P,
 768   gp_Vec& D1U, gp_Vec& D1V) const
 769 {
 770   if (equivSurf.IsNull()) {
 771     gp_Vec D2U,D2V,D2UV ;
 772     Handle(Geom_Surface) Basis = basisSurf;
 773
 774     // if Basis is Trimmed we take the basis's basis
 775     Handle(Geom_RectangularTrimmedSurface) RTS =
 776       Handle(Geom_RectangularTrimmedSurface)::DownCast(Basis);
 777     if (!RTS.IsNull())
 778       Basis = RTS->BasisSurface();
 779
 780     // BSpline case
 781     Handle(Geom_BSplineSurface) BSplS =
 782       Handle(Geom_BSplineSurface)::DownCast(Basis);
 783     if (!BSplS.IsNull()) {
 784       gp_Vec D3U,D3V,D3UUV,D3UVV;
 785       LocateSides(U,V,USide,VSide,BSplS,2,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 786       SetD1(U,V,P,D1U,D1V,D2U,D2V,D2UV);
 787       return;
 788     }
 789
 790     // Extrusion case
 791     Handle( Geom_SurfaceOfLinearExtrusion) SE =
 792       Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(Basis);
 793     if (!SE.IsNull()) {
 794       SE->LocalD2(U,V,USide,P,D1U,D1V,D2U,D2V,D2UV);
 795       SetD1(U,V,P,D1U,D1V,D2U,D2V,D2UV);
 796       return;
 797     }
 798
 799     // Revolution case
 800     Handle(Geom_SurfaceOfRevolution) SR =
 801       Handle(Geom_SurfaceOfRevolution)::DownCast(Basis);
 802     if (!SR.IsNull()) {
 803       SR->LocalD2(U,V,VSide,P,D1U,D1V,D2U,D2V,D2UV);
 804       SetD1(U,V,P,D1U,D1V,D2U,D2V,D2UV);
 805       return;
 806     }
 807
 808     // General cases
 809     basisSurf->D2(U, V, P, D1U, D1V, D2U, D2V, D2UV);
 810     SetD1(U,V,P,D1U,D1V,D2U,D2V,D2UV);
 811   }
 812   else
 813     equivSurf->D1(U,V,P,D1U,D1V);
 814 }
 815
 816 //=======================================================================
 817 //function : LocalD2
 818 //purpose  :
 819 //=======================================================================
 820
 821 void Geom_OffsetSurface::LocalD2 (const Standard_Real U, const Standard_Real V,
 822   const Standard_Integer USide, const Standard_Integer VSide,
 823   gp_Pnt& P,
 824   gp_Vec& D1U, gp_Vec& D1V,
 825   gp_Vec& D2U, gp_Vec& D2V, gp_Vec& D2UV) const
 826 {
 827   if (equivSurf.IsNull()) {
 828     Handle(Geom_Surface) Basis = basisSurf;
 829     gp_Vec D3U,D3V,D3UUV,D3UVV;
 830
 831     // if Basis is Trimmed we take the basis's basis
 832     Handle(Geom_RectangularTrimmedSurface) RTS =
 833       Handle(Geom_RectangularTrimmedSurface)::DownCast(Basis);
 834     if (!RTS.IsNull())
 835       Basis = RTS->BasisSurface();
 836
 837     // BSpline case
 838     Handle(Geom_BSplineSurface) BSplS =
 839       Handle(Geom_BSplineSurface)::DownCast(Basis);
 840     if (!BSplS.IsNull()) {
 841       LocateSides(U,V,USide,VSide,BSplS,3,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 842       SetD2(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 843       return;
 844     }
 845
 846     // Extrusion case
 847     Handle(Geom_SurfaceOfLinearExtrusion) SE =
 848       Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(Basis);
 849     if (!SE.IsNull()) {
 850       SE->LocalD3(U,V,USide,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 851       SetD2(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 852       return;
 853     }
 854
 855     // Revolution case
 856     Handle(Geom_SurfaceOfRevolution) SR =
 857       Handle(Geom_SurfaceOfRevolution)::DownCast(Basis);
 858     if (!SR.IsNull()) {
 859       SR->LocalD3(U,V,VSide,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 860       SetD2(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 861       return;
 862     }
 863
 864     // General cases
 865     basisSurf->D3(U, V, P, D1U, D1V, D2U, D2V, D2UV, D3U, D3V, D3UUV, D3UVV);
 866     SetD2(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 867   }
 868   else
 869     equivSurf->D2(U,V,P,D1U,D1V,D2U,D2V,D2UV);
 870 }
 871
 872 //=======================================================================
 873 //function : LocalD3
 874 //purpose  :
 875 //=======================================================================
 876
 877 void Geom_OffsetSurface::LocalD3 (const Standard_Real U, const Standard_Real V,
 878   const Standard_Integer USide, const Standard_Integer VSide,
 879   gp_Pnt& P,
 880   gp_Vec& D1U, gp_Vec& D1V,
 881   gp_Vec& D2U, gp_Vec& D2V, gp_Vec& D2UV,
 882   gp_Vec& D3U, gp_Vec& D3V, gp_Vec& D3UUV, gp_Vec& D3UVV) const
 883 {
 884   if (equivSurf.IsNull()) {
 885     Handle(Geom_Surface) Basis = basisSurf;
 886
 887     // if Basis is Trimmed we take the basis's basis
 888     Handle(Geom_RectangularTrimmedSurface) RTS =
 889       Handle(Geom_RectangularTrimmedSurface)::DownCast(Basis);
 890     if (!RTS.IsNull())
 891       Basis = RTS->BasisSurface();
 892
 893     // BSpline case
 894     Handle(Geom_BSplineSurface) BSplS =
 895       Handle(Geom_BSplineSurface)::DownCast(Basis);
 896     if (!BSplS.IsNull()) {
 897       LocateSides(U,V,USide,VSide,BSplS,3,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 898       SetD3(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 899       return;
 900     }
 901
 902     // Extrusion case
 903     Handle(Geom_SurfaceOfLinearExtrusion) SE =
 904       Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(Basis);
 905     if (!SE.IsNull()) {
 906       SE->LocalD3(U,V,USide,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 907       SetD3(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 908       return;
 909     }
 910
 911     // Revolution case
 912     Handle(Geom_SurfaceOfRevolution) SR =
 913       Handle(Geom_SurfaceOfRevolution)::DownCast(Basis);
 914     if (!SR.IsNull()) {
 915       SR->LocalD3(U,V,VSide,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 916       SetD3(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 917       return;
 918     }
 919
 920     // General cases
 921     basisSurf->D3(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 922     SetD3(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 923   }
 924   else
 925     equivSurf->D3(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 926 }
 927
 928 //=======================================================================
 929 //function : LocalDN
 930 //purpose  :
 931 //=======================================================================
 932
 933 gp_Vec Geom_OffsetSurface::LocalDN  (const Standard_Real U, const Standard_Real V,
 934   const Standard_Integer USide, const Standard_Integer VSide,
 935   const Standard_Integer Nu, const Standard_Integer Nv) const
 936 {
 937   if(equivSurf.IsNull()) {
 938
 939     gp_Vec D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV;
 940     gp_Pnt P;
 941     Handle(Geom_Surface) Basis = basisSurf;
 942
 943     // if Basis is Trimmed we make the basis`s basis
 944     Handle(Geom_RectangularTrimmedSurface) RTS =
 945       Handle(Geom_RectangularTrimmedSurface)::DownCast(Basis);
 946     if (!RTS.IsNull())
 947       Basis = RTS -> BasisSurface();
 948
 949     //BSpline case
 950     Handle(Geom_BSplineSurface) BSplS =
 951       Handle(Geom_BSplineSurface)::DownCast(Basis);
 952     if(!BSplS.IsNull()) {
 953       LocateSides(U,V,USide,VSide,BSplS,1,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV);
 954       return SetDN(U,V,Nu,Nv,D1U,D1V);
 955     }
 956
 957     //Extrusion case
 958     Handle( Geom_SurfaceOfLinearExtrusion) SE =
 959       Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(Basis);
 960     if(!SE.IsNull()) {
 961       SE->LocalD1(U,V,USide,P,D1U,D1V);
 962       return SetDN(U,V,Nu,Nv,D1U,D1V);
 963     }
 964
 965     //Revolution case
 966     Handle(Geom_SurfaceOfRevolution) SR =
 967       Handle(Geom_SurfaceOfRevolution)::DownCast(Basis);
 968     if(!SR.IsNull()) {
 969       SR->LocalDN(U,V,VSide,Nu,Nv);
 970       return SetDN(U,V,Nu,Nv,D1U,D1V);
 971     }
 972
 973     //General cases
 974     basisSurf->DN(U,V,Nu,Nv);
 975     return SetDN(U,V,Nu,Nv,D1U,D1V);
 976   }
 977   else
 978     return equivSurf->DN(U,V,Nu,Nv);
 979 }
 980
 981 //=======================================================================
 982 //function : LocateSides
 983 //purpose  :
 984 //=======================================================================
 985
 986 void Geom_OffsetSurface::LocateSides(const Standard_Real U, const Standard_Real V,
 987   const Standard_Integer USide, const Standard_Integer VSide,
 988   const Handle(Geom_BSplineSurface)& BSplS,
 989   const Standard_Integer NDir,
 990   gp_Pnt& P,
 991   gp_Vec& D1U,gp_Vec& D1V,
 992   gp_Vec& D2U,gp_Vec& D2V,gp_Vec& D2UV,
 993   gp_Vec& D3U,gp_Vec& D3V,gp_Vec& D3UUV,gp_Vec& D3UVV) const
 994 {
 995   Standard_Boolean UIsKnot=Standard_False, VIsKnot=Standard_False;
 996   Standard_Integer Ideb, Ifin, IVdeb, IVfin;
 997   const Standard_Real ParTol=Precision::PConfusion()/10.;
 998   BSplS->Geom_BSplineSurface::LocateU(U,ParTol,Ideb, Ifin, Standard_False);
 999   BSplS->Geom_BSplineSurface::LocateV(V,ParTol,IVdeb,IVfin,Standard_False);
1000
1001   if(Ideb == Ifin ) { //knot
1002
1003     if(USide == 1) {Ifin++; UIsKnot=Standard_True;}
1004     else if(USide == -1) {Ideb--; UIsKnot=Standard_True;}
1005     else {Ideb--; Ifin++;} //USide == 0
1006
1007   }
1008
1009   if(Ideb < BSplS->FirstUKnotIndex()) {Ideb = BSplS->FirstUKnotIndex(); Ifin = Ideb + 1;}
1010
1011   if(Ifin > BSplS->LastUKnotIndex()) {Ifin = BSplS->LastUKnotIndex(); Ideb = Ifin - 1;}
1012
1013   if(IVdeb == IVfin ) { //knot
1014
1015     if(VSide == 1) {IVfin++; VIsKnot=Standard_True;}
1016     else if(VSide == -1) {IVdeb--; VIsKnot=Standard_True;}
1017     else {IVdeb--; IVfin++;} //VSide == 0
1018
1019   }
1020
1021   if(IVdeb < BSplS->FirstVKnotIndex()) {IVdeb = BSplS->FirstVKnotIndex(); IVfin = IVdeb + 1;}
1022
1023   if(IVfin > BSplS->LastVKnotIndex()) {IVfin = BSplS->LastVKnotIndex(); IVdeb = IVfin - 1;}
1024
1025   if((UIsKnot)||(VIsKnot))
1026     switch(NDir)
1027     {
1028       case 0 : BSplS->Geom_BSplineSurface::LocalD0(U,V,Ideb,Ifin,IVdeb,IVfin,P); break;
1029       case 1 : BSplS->Geom_BSplineSurface::LocalD1(U,V,Ideb,Ifin,IVdeb,IVfin,P,D1U,D1V); break;
1030       case 2 : BSplS->Geom_BSplineSurface::LocalD2(U,V,Ideb,Ifin,IVdeb,IVfin,P,D1U,D1V,D2U,D2V,D2UV); break;
1031       case 3 : BSplS->Geom_BSplineSurface::LocalD3(U,V,Ideb,Ifin,IVdeb,IVfin,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV); break;
1032     }
1033   else
1034     switch(NDir)
1035     {
1036       case 0 : basisSurf->D0(U,V,P); break;
1037       case 1 : basisSurf->D1(U,V,P,D1U,D1V); break;
1038       case 2 : basisSurf->D2(U,V,P,D1U,D1V,D2U,D2V,D2UV); break;
1039       case 3 : basisSurf->D3(U,V,P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV); break;
1040     }
1041 }
1042
1043 ////*************************************************
1044 ////
1045 ////   EVALUATOR FOR THE ISO-CURVE APPROXIMATION
1046 ////
1047 ////*************************************************
1048
1049 class Geom_OffsetSurface_UIsoEvaluator : public AdvApprox_EvaluatorFunction
1050 {
1051 public:
1052   Geom_OffsetSurface_UIsoEvaluator (const Handle(Geom_Surface)& theSurface, const Standard_Real theU)
1053     : CurrentSurface(theSurface), IsoPar(theU) {}
1054
1055   virtual void Evaluate (Standard_Integer *Dimension,
1056     Standard_Real     StartEnd[2],
1057     Standard_Real    *Parameter,
1058     Standard_Integer *DerivativeRequest,
1059     Standard_Real    *Result, // [Dimension]
1060     Standard_Integer *ErrorCode);
1061
1062 private:
1063   Handle(Geom_Surface) CurrentSurface;
1064   Standard_Real IsoPar;
1065 };
1066
1067 void Geom_OffsetSurface_UIsoEvaluator::Evaluate(Standard_Integer *,/*Dimension*/
1068   Standard_Real     /*StartEnd*/[2],
1069   Standard_Real    *Parameter,
1070   Standard_Integer *DerivativeRequest,
1071   Standard_Real    *Result,
1072   Standard_Integer *ReturnCode)
1073 {
1074   gp_Pnt P;
1075   if (*DerivativeRequest == 0) {
1076     P = CurrentSurface->Value(IsoPar,*Parameter);
1077     Result[0] = P.X();
1078     Result[1] = P.Y();
1079     Result[2] = P.Z();
1080   }
1081   else {
1082     gp_Vec DU,DV;
1083     CurrentSurface->D1(IsoPar,*Parameter,P,DU,DV);
1084     Result[0] = DV.X();
1085     Result[1] = DV.Y();
1086     Result[2] = DV.Z();
1087   }
1088   *ReturnCode = 0;
1089 }
1090
1091 class Geom_OffsetSurface_VIsoEvaluator : public AdvApprox_EvaluatorFunction
1092 {
1093 public:
1094   Geom_OffsetSurface_VIsoEvaluator (const Handle(Geom_Surface)& theSurface, const Standard_Real theV)
1095     : CurrentSurface(theSurface), IsoPar(theV) {}
1096
1097   virtual void Evaluate (Standard_Integer *Dimension,
1098     Standard_Real     StartEnd[2],
1099     Standard_Real    *Parameter,
1100     Standard_Integer *DerivativeRequest,
1101     Standard_Real    *Result, // [Dimension]
1102     Standard_Integer *ErrorCode);
1103
1104 private:
1105   Handle(Geom_Surface) CurrentSurface;
1106   Standard_Real IsoPar;
1107 };
1108
1109 void Geom_OffsetSurface_VIsoEvaluator::Evaluate(Standard_Integer *,/*Dimension*/
1110   Standard_Real     /*StartEnd*/[2],
1111   Standard_Real    *Parameter,
1112   Standard_Integer *DerivativeRequest,
1113   Standard_Real    *Result,
1114   Standard_Integer *ReturnCode)
1115 {
1116   gp_Pnt P;
1117   if (*DerivativeRequest == 0) {
1118     P = CurrentSurface->Value(*Parameter,IsoPar);
1119     Result[0] = P.X();
1120     Result[1] = P.Y();
1121     Result[2] = P.Z();
1122   }
1123   else {
1124     gp_Vec DU,DV;
1125     CurrentSurface->D1(*Parameter,IsoPar,P,DU,DV);
1126     Result[0] = DU.X();
1127     Result[1] = DU.Y();
1128     Result[2] = DU.Z();
1129   }
1130   *ReturnCode = 0;
1131 }
1132
1133 //=======================================================================
1134 //function : UIso
1135 //purpose  : The Uiso or the VIso of an OffsetSurface can't be clearly
1136 //           exprimed as a curve from Geom. So, to extract the U or VIso
1137 //           an Approximation is needed. This approx always will return a
1138 //           BSplineCurve from Geom.
1139 //=======================================================================
1140
1141 Handle(Geom_Curve) Geom_OffsetSurface::UIso (const Standard_Real UU) const
1142 {
1143   if (equivSurf.IsNull()) {
1144     const Standard_Integer Num1 = 0, Num2 = 0, Num3 = 1;
1145     Handle(TColStd_HArray1OfReal) T1, T2, T3 = new TColStd_HArray1OfReal(1,Num3);
1146     T3->Init(Precision::Approximation());
1147     Standard_Real U1,U2,V1,V2;
1148     Bounds(U1,U2,V1,V2);
1149     const GeomAbs_Shape Cont = GeomAbs_C1;
1150     const Standard_Integer MaxSeg = 100, MaxDeg = 14;
1151
1152     Handle(Geom_OffsetSurface) me (this);
1153     Geom_OffsetSurface_UIsoEvaluator ev (me, UU);
1154     AdvApprox_ApproxAFunction Approx(Num1, Num2, Num3, T1, T2, T3,
1155       V1, V2, Cont, MaxDeg, MaxSeg, ev);
1156
1157     Standard_ConstructionError_Raise_if (!Approx.IsDone(), " Geom_OffsetSurface : UIso");
1158
1159     const Standard_Integer NbPoles = Approx.NbPoles();
1160
1161     TColgp_Array1OfPnt      Poles( 1, NbPoles);
1162     TColStd_Array1OfReal    Knots( 1, Approx.NbKnots());
1163     TColStd_Array1OfInteger Mults( 1, Approx.NbKnots());
1164
1165     Approx.Poles(1, Poles);
1166     Knots = Approx.Knots()->Array1();
1167     Mults = Approx.Multiplicities()->Array1();
1168
1169     Handle(Geom_BSplineCurve) C =
1170       new Geom_BSplineCurve( Poles, Knots, Mults, Approx.Degree());
1171     return C;
1172   }
1173   else
1174     return equivSurf->UIso(UU);
1175 }
1176
1177 //=======================================================================
1178 //function : Value
1179 //purpose  :
1180 //=======================================================================
1181
1182 void Geom_OffsetSurface::Value
1183   (const Standard_Real U,  const Standard_Real V,
1184    gp_Pnt& P,        gp_Pnt& ,
1185    gp_Vec& D1Ubasis, gp_Vec& D1Vbasis) const
1186 {
1187   if (myBasisSurfContinuity == GeomAbs_C0)
1188     Geom_UndefinedValue::Raise();
1189
1190   SetD0(U,V,P,D1Ubasis,D1Vbasis);
1191 }
1192
1193 //=======================================================================
1194 //function : VIso
1195 //purpose  :
1196 //=======================================================================
1197
1198 Handle(Geom_Curve) Geom_OffsetSurface::VIso (const Standard_Real VV) const
1199 {
1200   if (equivSurf.IsNull()) {
1201     const Standard_Integer Num1 = 0, Num2 = 0, Num3 = 1;
1202     Handle(TColStd_HArray1OfReal) T1, T2, T3 = new TColStd_HArray1OfReal(1,Num3);
1203     T3->Init(Precision::Approximation());
1204     Standard_Real U1,U2,V1,V2;
1205     Bounds(U1,U2,V1,V2);
1206     const GeomAbs_Shape Cont = GeomAbs_C1;
1207     const Standard_Integer MaxSeg = 100, MaxDeg = 14;
1208
1209     Handle(Geom_OffsetSurface) me (this);
1210     Geom_OffsetSurface_VIsoEvaluator ev (me, VV);
1211     AdvApprox_ApproxAFunction Approx (Num1, Num2, Num3, T1, T2, T3,
1212       U1, U2, Cont, MaxDeg, MaxSeg, ev);
1213
1214     Standard_ConstructionError_Raise_if (!Approx.IsDone(), " Geom_OffsetSurface : VIso");
1215
1216     TColgp_Array1OfPnt      Poles( 1, Approx.NbPoles());
1217     TColStd_Array1OfReal    Knots( 1, Approx.NbKnots());
1218     TColStd_Array1OfInteger Mults( 1, Approx.NbKnots());
1219
1220     Approx.Poles(1, Poles);
1221     Knots = Approx.Knots()->Array1();
1222     Mults = Approx.Multiplicities()->Array1();
1223
1224     Handle(Geom_BSplineCurve) C =
1225       new Geom_BSplineCurve( Poles, Knots, Mults, Approx.Degree());
1226     return C;
1227   }
1228   else
1229     return equivSurf->VIso(VV);
1230 }
1231
1232 //=======================================================================
1233 //function : IsCNu
1234 //purpose  :
1235 //=======================================================================
1236
1237 Standard_Boolean Geom_OffsetSurface::IsCNu (const Standard_Integer N) const
1238 {
1239   Standard_RangeError_Raise_if (N < 0, " ");
1240   return basisSurf->IsCNu (N+1);
1241 }
1242
1243 //=======================================================================
1244 //function : IsCNv
1245 //purpose  :
1246 //=======================================================================
1247
1248 Standard_Boolean Geom_OffsetSurface::IsCNv (const Standard_Integer N) const
1249 {
1250   Standard_RangeError_Raise_if (N < 0, " ");
1251   return basisSurf->IsCNv (N+1);
1252 }
1253
1254 //=======================================================================
1255 //function : IsUPeriodic
1256 //purpose  :
1257 //=======================================================================
1258
1259 Standard_Boolean Geom_OffsetSurface::IsUPeriodic () const
1260 {
1261   return basisSurf->IsUPeriodic();
1262 }
1263
1264 //=======================================================================
1265 //function : UPeriod
1266 //purpose  :
1267 //=======================================================================
1268
1269 Standard_Real Geom_OffsetSurface::UPeriod() const
1270 {
1271   return basisSurf->UPeriod();
1272 }
1273
1274 //=======================================================================
1275 //function : IsVPeriodic
1276 //purpose  :
1277 //=======================================================================
1278
1279 Standard_Boolean Geom_OffsetSurface::IsVPeriodic () const
1280 {
1281   return basisSurf->IsVPeriodic();
1282 }
1283
1284 //=======================================================================
1285 //function : VPeriod
1286 //purpose  :
1287 //=======================================================================
1288
1289 Standard_Real Geom_OffsetSurface::VPeriod() const
1290 {
1291   return basisSurf->VPeriod();
1292 }
1293
1294 //=======================================================================
1295 //function : IsUClosed
1296 //purpose  :
1297 //=======================================================================
1298
1299 Standard_Boolean Geom_OffsetSurface::IsUClosed () const
1300 {
1301   Standard_Boolean UClosed;
1302   Handle(Geom_Surface) SBasis = BasisSurface();
1303
1304   if (SBasis->IsKind (STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
1305     Handle(Geom_RectangularTrimmedSurface) St =
1306       Handle(Geom_RectangularTrimmedSurface)::DownCast(SBasis);
1307
1308     Handle(Geom_Surface) S = Handle(Geom_Surface)::DownCast(St->BasisSurface());
1309     if (S->IsKind (STANDARD_TYPE(Geom_ElementarySurface))) {
1310       UClosed = SBasis->IsUClosed();
1311     }
1312     else if (S->IsKind (STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion))) {
1313       Handle(Geom_SurfaceOfLinearExtrusion) Extru =
1314         Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(S);
1315
1316       Handle(Geom_Curve) C = Extru->BasisCurve();
1317       if (C->IsKind (STANDARD_TYPE(Geom_Circle)) || C->IsKind (STANDARD_TYPE(Geom_Ellipse))) {
1318         UClosed = SBasis->IsUClosed();
1319       }
1320       else { UClosed = Standard_False; }
1321     }
1322     else if (S->IsKind (STANDARD_TYPE(Geom_SurfaceOfRevolution))) {
1323       UClosed = SBasis->IsUClosed();
1324     }
1325     else { UClosed = Standard_False; }
1326   }
1327   else {
1328     if (SBasis->IsKind (STANDARD_TYPE(Geom_ElementarySurface))) {
1329       UClosed = SBasis->IsUClosed();
1330     }
1331     else if (SBasis->IsKind (STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion))) {
1332       Handle(Geom_SurfaceOfLinearExtrusion) Extru =
1333         Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(SBasis);
1334
1335       Handle(Geom_Curve) C = Extru->BasisCurve();
1336       UClosed = (C->IsKind(STANDARD_TYPE(Geom_Circle)) || C->IsKind(STANDARD_TYPE(Geom_Ellipse)));
1337     }
1338     else if (SBasis->IsKind (STANDARD_TYPE(Geom_SurfaceOfRevolution))) {
1339       UClosed = Standard_True;
1340     }
1341     else { UClosed = Standard_False; }
1342   }
1343   return UClosed;
1344 }
1345
1346 //=======================================================================
1347 //function : IsVClosed
1348 //purpose  :
1349 //=======================================================================
1350
1351 Standard_Boolean Geom_OffsetSurface::IsVClosed () const
1352 {
1353   Standard_Boolean VClosed;
1354   Handle(Geom_Surface) SBasis = BasisSurface();
1355
1356   if (SBasis->IsKind (STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
1357     Handle(Geom_RectangularTrimmedSurface) St =
1358       Handle(Geom_RectangularTrimmedSurface)::DownCast(SBasis);
1359
1360     Handle(Geom_Surface) S = Handle(Geom_Surface)::DownCast(St->BasisSurface());
1361     if (S->IsKind (STANDARD_TYPE(Geom_ElementarySurface))) {
1362       VClosed = SBasis->IsVClosed();
1363     }
1364     else { VClosed = Standard_False; }
1365   }
1366   else {
1367     if (SBasis->IsKind (STANDARD_TYPE(Geom_ElementarySurface))) {
1368       VClosed = SBasis->IsVClosed();
1369     }
1370     else { VClosed = Standard_False; }
1371   }
1372   return VClosed;
1373 }
1374
1375 //=======================================================================
1376 //function : Transform
1377 //purpose  :
1378 //=======================================================================
1379
1380 void Geom_OffsetSurface::Transform (const gp_Trsf& T)
1381 {
1382   basisSurf->Transform (T);
1383   offsetValue *= T.ScaleFactor();
1384   equivSurf.Nullify();
1385 }
1386
1387 //=======================================================================
1388 //function : TransformParameters
1389 //purpose  :
1390 //=======================================================================
1391
1392 void Geom_OffsetSurface::TransformParameters(Standard_Real& U, Standard_Real& V,
1393   const gp_Trsf& T) const
1394 {
1395   basisSurf->TransformParameters(U,V,T);
1396   if(!equivSurf.IsNull()) equivSurf->TransformParameters(U,V,T);
1397 }
1398
1399 //=======================================================================
1400 //function : ParametricTransformation
1401 //purpose  :
1402 //=======================================================================
1403
1404 gp_GTrsf2d Geom_OffsetSurface::ParametricTransformation (const gp_Trsf& T) const
1405 {
1406   return basisSurf->ParametricTransformation(T);
1407 }
1408
1409 //=======================================================================
1410 //function : Surface
1411 //purpose  : Trouve si elle existe, une surface non offset, equivalente
1412 //           a l'offset surface.
1413 //=======================================================================
1414
1415 Handle(Geom_Surface) Geom_OffsetSurface::Surface() const
1416 {
1417   if (offsetValue == 0.0) return  basisSurf; // Cas direct
1418
1419   Standard_Real Tol = Precision::Confusion();
1420   Handle(Geom_Surface) Result, Base;
1421   Result.Nullify();
1422   Handle(Standard_Type) TheType = basisSurf->DynamicType();
1423   Standard_Boolean IsTrimmed;
1424   Standard_Real U1 = 0., V1 = 0., U2 = 0., V2 = 0.;
1425
1426   // Preambule pour les surface trimmes
1427   if (TheType == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
1428     Handle(Geom_RectangularTrimmedSurface) S =
1429       Handle(Geom_RectangularTrimmedSurface)::DownCast(basisSurf);
1430     Base = S->BasisSurface();
1431     TheType = Base->DynamicType();
1432     S->Bounds(U1,U2,V1,V2);
1433     IsTrimmed = Standard_True;
1434   }
1435   else {
1436     IsTrimmed = Standard_False;
1437     Base = basisSurf;
1438   }
1439
1440   // Traite les surfaces cannonique
1441   if (TheType == STANDARD_TYPE(Geom_Plane))
1442   {
1443     Handle(Geom_Plane) P =
1444       Handle(Geom_Plane)::DownCast(Base);
1445     gp_Vec T = P->Position().XDirection()^P->Position().YDirection();
1446     T *= offsetValue;
1447     Result = Handle(Geom_Plane)::DownCast(P->Translated(T));
1448   }
1449   else if (TheType == STANDARD_TYPE(Geom_CylindricalSurface))
1450   {
1451     Handle(Geom_CylindricalSurface) C =
1452       Handle(Geom_CylindricalSurface)::DownCast(Base);
1453     Standard_Real Radius = C->Radius();
1454     gp_Ax3 Axis = C->Position();
1455     if (Axis.Direct())
1456       Radius += offsetValue;
1457     else
1458       Radius -= offsetValue;
1459     if ( Radius >= Tol ) {
1460       Result = new Geom_CylindricalSurface( Axis, Radius);
1461     }
1462     else if ( Radius <= -Tol ){
1463       Axis.Rotate(gp_Ax1(Axis.Location(),Axis.Direction()),M_PI);
1464       Result = new Geom_CylindricalSurface( Axis, Abs(Radius));
1465       Result->UReverse();
1466     }
1467     else
1468     {
1469       // surface degeneree
1470     }
1471   }
1472   else if (TheType == STANDARD_TYPE(Geom_ConicalSurface))
1473   {
1474     Handle(Geom_ConicalSurface) C =
1475       Handle(Geom_ConicalSurface)::DownCast(Base);
1476     gp_Ax3 anAxis = C->Position();
1477     Standard_Boolean isDirect = anAxis.Direct();
1478     Standard_Real anAlpha = C->SemiAngle();
1479     Standard_Real aRadius;
1480     if (isDirect)
1481     {
1482       aRadius = C->RefRadius() + offsetValue * Cos (anAlpha);
1483     }
1484     else
1485     {
1486       aRadius = C->RefRadius() - offsetValue * Cos (anAlpha);
1487     }
1488     if (aRadius >= 0.)
1489     {
1490       gp_Vec aZ (anAxis.Direction());
1491       if (isDirect)
1492       {
1493         aZ *= -offsetValue * Sin (anAlpha);
1494       }
1495       else
1496       {
1497         aZ *=  offsetValue * Sin (anAlpha);
1498       }
1499       anAxis.Translate (aZ);
1500       Result = new Geom_ConicalSurface (anAxis, anAlpha, aRadius);
1501     }
1502     else
1503     {
1504       // surface degeneree
1505     }
1506   }
1507   else if (TheType == STANDARD_TYPE(Geom_SphericalSurface)) {
1508     Handle(Geom_SphericalSurface) S =
1509       Handle(Geom_SphericalSurface)::DownCast(Base);
1510     Standard_Real Radius = S->Radius();
1511     gp_Ax3 Axis = S->Position();
1512     if (Axis.Direct())
1513       Radius += offsetValue;
1514     else
1515       Radius -= offsetValue;
1516     if ( Radius >= Tol) {
1517       Result = new Geom_SphericalSurface(Axis, Radius);
1518     }
1519     else if ( Radius <= -Tol ) {
1520       Axis.Rotate(gp_Ax1(Axis.Location(),Axis.Direction()),M_PI);
1521       Axis.ZReverse();
1522       Result = new Geom_SphericalSurface(Axis, -Radius);
1523       Result->UReverse();
1524     }
1525     else {
1526       //      surface degeneree
1527     }
1528   }
1529   else if (TheType == STANDARD_TYPE(Geom_ToroidalSurface))
1530
1531   {
1532     Handle(Geom_ToroidalSurface)
1533       S = Handle(Geom_ToroidalSurface)::DownCast(Base);
1534     Standard_Real MajorRadius = S->MajorRadius();
1535     Standard_Real MinorRadius = S->MinorRadius();
1536     gp_Ax3 Axis = S->Position();
1537     if (MinorRadius <= MajorRadius)
1538     {
1539       if (Axis.Direct())
1540         MinorRadius += offsetValue;
1541       else
1542         MinorRadius -= offsetValue;
1543       if (MinorRadius >= Tol)
1544         Result = new Geom_ToroidalSurface(Axis,MajorRadius,MinorRadius);
1545       //      else if (MinorRadius <= -Tol)
1546       //        Result->UReverse();
1547       else
1548       {
1549         //      surface degeneree
1550       }
1551     }
1552   }
1553
1554   // S'il le faut on trimme le resultat
1555   if (IsTrimmed && !Result.IsNull()) {
1556     Base = Result;
1557     Result = new Geom_RectangularTrimmedSurface (Base, U1, U2, V1,V2);
1558   }
1559
1560   return Result;
1561 }
1562
1563 //=======================================================================
1564 //function : UOsculatingSurface
1565 //purpose  :
1566 //=======================================================================
1567
1568 Standard_Boolean Geom_OffsetSurface::UOsculatingSurface(const Standard_Real U, const Standard_Real V,
1569   Standard_Boolean& t, Handle(Geom_BSplineSurface)& L) const
1570 {
1571   return myOscSurf.UOscSurf(U,V,t,L);
1572 }
1573
1574 //=======================================================================
1575 //function : VOsculatingSurface
1576 //purpose  :
1577 //=======================================================================
1578
1579 Standard_Boolean Geom_OffsetSurface::VOsculatingSurface(const Standard_Real U, const Standard_Real V,
1580   Standard_Boolean& t, Handle(Geom_BSplineSurface)& L) const
1581 {
1582   return myOscSurf.VOscSurf(U,V,t,L);
1583 }
1584
1585 //=======================================================================
1586 //function :
1587 //purpose  : private
1588 //=======================================================================
1589
1590 void Geom_OffsetSurface::SetD0(const Standard_Real U, const Standard_Real V,
1591                                gp_Pnt& P,
1592                                const gp_Vec& D1U, const gp_Vec& D1V) const // First Derivative
1593 {
1594   const Standard_Real MagTol=0.000000001;
1595
1596   gp_Vec aNorm = D1U.Crossed(D1V);
1597   if (aNorm.SquareMagnitude() > MagTol * MagTol)
1598   {
1599     // Non singular case. Simple computations.
1600     aNorm.Normalize();
1601     P.SetXYZ(P.XYZ() + offsetValue * aNorm.XYZ());
1602   }
1603   else
1604   {
1605     Standard_Boolean AlongU = Standard_False,
1606                      AlongV = Standard_False;
1607     Handle(Geom_BSplineSurface) L;
1608     Standard_Boolean IsOpposite=Standard_False;
1609     AlongU = UOsculatingSurface(U,V,IsOpposite,L);
1610     AlongV = VOsculatingSurface(U,V,IsOpposite,L);
1611     gp_Dir Normal;
1612     CSLib_NormalStatus NStatus;
1613     Standard_Integer MaxOrder=3;
1614     TColgp_Array2OfVec DerNUV(0,MaxOrder,0,MaxOrder);
1615     TColgp_Array2OfVec DerSurf(0,MaxOrder+1,0,MaxOrder+1);
1616     Standard_Integer OrderU,OrderV;
1617     Standard_Real Umin,Umax,Vmin,Vmax;
1618     Bounds(Umin,Umax,Vmin,Vmax);
1619     DerSurf.SetValue(1, 0, D1U);
1620     DerSurf.SetValue(0, 1, D1V);
1621     derivatives(MaxOrder,1,U,V,basisSurf,0,0,AlongU,AlongV,L,DerNUV,DerSurf);
1622
1623     Standard_Real signe = 1.0;
1624     if ((AlongV || AlongU) && IsOpposite)
1625       signe = -1;
1626
1627     CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,NStatus,Normal,OrderU,OrderV);
1628     if (NStatus == CSLib_Defined)
1629       P.SetXYZ(P.XYZ() + offsetValue * signe * Normal.XYZ());
1630     else
1631     {
1632       if (NStatus == CSLib_InfinityOfSolutions &&
1633           D1U.SquareMagnitude() + D1V.SquareMagnitude() > MagTol * MagTol)
1634       {
1635         // Use non-null derivative as normal direction in degenerated case.
1636         aNorm = D1U.SquareMagnitude() > MagTol * MagTol ? D1U : D1V;
1637         aNorm.Normalize();
1638
1639         P.SetXYZ(P.XYZ() + offsetValue * signe * aNorm.XYZ());
1640         return;
1641       }
1642       Geom_UndefinedValue::Raise();
1643     }
1644   }
1645 }
1646
1647 //=======================================================================
1648 //function :
1649 //purpose  : private
1650 //=======================================================================/
1651
1652 void Geom_OffsetSurface::SetD1(const Standard_Real U, const Standard_Real V,
1653                                gp_Pnt& P,
1654                                gp_Vec& D1U, gp_Vec& D1V, // First derivative
1655                                const gp_Vec& D2UU, const gp_Vec& D2VV, const gp_Vec& D2UV) const // Second derivative
1656 {
1657   const Standard_Real MagTol=0.000000001;
1658
1659   // Check offset side.
1660   Handle(Geom_BSplineSurface) L;
1661   Standard_Boolean AlongU = Standard_False,
1662                    AlongV = Standard_False;
1663   Standard_Boolean IsOpposite=Standard_False;
1664   AlongU = UOsculatingSurface(U,V,IsOpposite,L);
1665   AlongV = VOsculatingSurface(U,V,IsOpposite,L);
1666   Standard_Real signe = 1.0;
1667   if ((AlongV || AlongU) && IsOpposite)
1668     signe = -1.0;
1669
1670   Standard_Integer MaxOrder = 0;
1671   gp_Vec aNorm = D1U.Crossed(D1V);
1672   if (aNorm.SquareMagnitude() > MagTol * MagTol)
1673   {
1674     MaxOrder=0;
1675
1676     if (!AlongV && !AlongU)
1677     {
1678       // AlongU or AlongV leads to more complex D1 computations,
1679       // try to compute D0 and D1 much simpler.
1680       aNorm.Normalize();
1681       P.SetXYZ(P.XYZ() + offsetValue * signe * aNorm.XYZ());
1682
1683       gp_Vec aN0(aNorm.XYZ()), aN1U, aN1V;
1684       Standard_Real aScale = (D1U^D1V).Dot(aN0);
1685       aN1U.SetX(D2UU.Y() * D1V.Z() + D1U.Y() * D2UV.Z()
1686               - D2UU.Z() * D1V.Y() - D1U.Z() * D2UV.Y());
1687       aN1U.SetY((D2UU.X() * D1V.Z() + D1U.X() * D2UV.Z()
1688                - D2UU.Z() * D1V.X() - D1U.Z() * D2UV.X() ) * -1.0);
1689       aN1U.SetZ(D2UU.X() * D1V.Y() + D1U.X() * D2UV.Y()
1690               - D2UU.Y() * D1V.X() - D1U.Y() * D2UV.X());
1691       Standard_Real aScaleU = aN1U.Dot(aN0);
1692       aN1U.Subtract(aScaleU * aN0);
1693       aN1U /= aScale;
1694
1695       aN1V.SetX(D2UV.Y() * D1V.Z() + D2VV.Z() * D1U.Y()
1696               - D2UV.Z() * D1V.Y() - D2VV.Y() * D1U.Z());
1697       aN1V.SetY((D2UV.X() * D1V.Z() + D2VV.Z() * D1U.X()
1698                - D2UV.Z() * D1V.X() - D2VV.X() * D1U.Z()) * -1.0);
1699       aN1V.SetZ(D2UV.X() * D1V.Y() + D2VV.Y() * D1U.X()
1700               - D2UV.Y() * D1V.X() - D2VV.X() * D1U.Y());
1701       Standard_Real aScaleV = aN1V.Dot(aN0);
1702       aN1V.Subtract(aScaleV * aN0);
1703       aN1V /= aScale;
1704
1705       D1U += offsetValue * signe * aN1U;
1706       D1V += offsetValue * signe * aN1V;
1707
1708       return;
1709     }
1710   }
1711   else
1712     MaxOrder=3;
1713
1714   Standard_Integer OrderU,OrderV;
1715   TColgp_Array2OfVec DerNUV(0,MaxOrder+1,0,MaxOrder+1);
1716   TColgp_Array2OfVec DerSurf(0,MaxOrder+2,0,MaxOrder+2);
1717   Standard_Real Umin,Umax,Vmin,Vmax;
1718   Bounds(Umin,Umax,Vmin,Vmax);
1719   DerSurf.SetValue(1, 0, D1U);
1720   DerSurf.SetValue(0, 1, D1V);
1721   DerSurf.SetValue(1, 1, D2UV);
1722   DerSurf.SetValue(2, 0, D2UU);
1723   DerSurf.SetValue(0, 2, D2VV);
1724   derivatives(MaxOrder,2,U,V,basisSurf,1,1,AlongU,AlongV,L,DerNUV,DerSurf);
1725
1726   gp_Dir Normal;
1727   CSLib_NormalStatus NStatus;
1728   CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,NStatus,Normal,OrderU,OrderV);
1729   if (NStatus != CSLib_Defined) Geom_UndefinedValue::Raise();
1730
1731   P.SetXYZ(P.XYZ() + offsetValue * signe * Normal.XYZ());
1732
1733   D1U = DerSurf(1,0)
1734     + offsetValue * signe * CSLib::DNNormal(1,0,DerNUV,OrderU,OrderV);
1735   D1V = DerSurf(0,1)
1736     + offsetValue * signe * CSLib::DNNormal(0,1,DerNUV,OrderU,OrderV);
1737 }
1738
1739 //=======================================================================
1740 //function :
1741 //purpose  : private
1742 //=======================================================================/
1743
1744 void Geom_OffsetSurface::SetD2(const Standard_Real U, const Standard_Real V,
1745   gp_Pnt& P,
1746   gp_Vec& D1U, gp_Vec& D1V,
1747   gp_Vec& D2U, gp_Vec& D2V, gp_Vec& D2UV,
1748   const gp_Vec& d3u, const gp_Vec& d3v,
1749   const gp_Vec& d3uuv, const gp_Vec& d3uvv) const
1750 {
1751   const Standard_Real MagTol=0.000000001;
1752
1753   gp_Dir Normal;
1754   CSLib_NormalStatus NStatus;
1755   CSLib::Normal (D1U, D1V, MagTol, NStatus, Normal);
1756
1757   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined)? 0 : 3;
1758   Standard_Integer OrderU,OrderV;
1759   TColgp_Array2OfVec DerNUV(0,MaxOrder+2,0,MaxOrder+2);
1760   TColgp_Array2OfVec DerSurf(0,MaxOrder+3,0,MaxOrder+3);
1761   Standard_Real Umin,Umax,Vmin,Vmax;
1762   Bounds(Umin,Umax,Vmin,Vmax);
1763   DerSurf.SetValue(1, 0, D1U);
1764   DerSurf.SetValue(0, 1, D1V);
1765   DerSurf.SetValue(1, 1, D2UV);
1766   DerSurf.SetValue(2, 0, D2U);
1767   DerSurf.SetValue(0, 2, D2V);
1768   DerSurf.SetValue(3, 0, d3u);
1769   DerSurf.SetValue(2, 1, d3uuv);
1770   DerSurf.SetValue(1, 2, d3uvv);
1771   DerSurf.SetValue(0, 3, d3v);
1772   //*********************
1773
1774   Handle(Geom_BSplineSurface) L;
1775   Standard_Boolean IsOpposite=Standard_False;
1776   const Standard_Boolean AlongU = UOsculatingSurface(U,V,IsOpposite,L);
1777   const Standard_Boolean AlongV = VOsculatingSurface(U,V,IsOpposite,L);
1778   const Standard_Real signe = ((AlongV || AlongU) && IsOpposite)? -1. : 1.;
1779   derivatives(MaxOrder,3,U,V,basisSurf,2,2,AlongU,AlongV,L,DerNUV,DerSurf);
1780
1781   CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,NStatus,Normal,OrderU,OrderV);
1782   if (NStatus != CSLib_Defined) Geom_UndefinedValue::Raise();
1783
1784   P.SetXYZ(P.XYZ() + offsetValue * signe * Normal.XYZ());
1785
1786   D1U = DerSurf(1,0)
1787     + offsetValue * signe * CSLib::DNNormal(1,0,DerNUV,OrderU,OrderV);
1788   D1V = DerSurf(0,1)
1789     + offsetValue * signe * CSLib::DNNormal(0,1,DerNUV,OrderU,OrderV);
1790
1791   D2U = basisSurf->DN(U,V,2,0)
1792     + signe * offsetValue * CSLib::DNNormal(2,0,DerNUV,OrderU,OrderV);
1793   D2V = basisSurf->DN(U,V,0,2)
1794     + signe * offsetValue * CSLib::DNNormal(0,2,DerNUV,OrderU,OrderV);
1795   D2UV = basisSurf->DN(U,V,1,1)
1796     + signe * offsetValue * CSLib::DNNormal(1,1,DerNUV,OrderU,OrderV);
1797 }
1798
1799 //=======================================================================
1800 //function :
1801 //purpose  : private
1802 //=======================================================================/
1803
1804 void Geom_OffsetSurface::SetD3(const Standard_Real U, const Standard_Real V,
1805   gp_Pnt& P,
1806   gp_Vec& D1U, gp_Vec& D1V,
1807   gp_Vec& D2U, gp_Vec& D2V, gp_Vec& D2UV,
1808   gp_Vec& D3U, gp_Vec& D3V, gp_Vec& D3UUV, gp_Vec& D3UVV) const
1809 {
1810   const Standard_Real MagTol=0.000000001;
1811
1812   gp_Dir Normal;
1813   CSLib_NormalStatus NStatus;
1814   CSLib::Normal (D1U, D1V, MagTol, NStatus, Normal);
1815   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined)? 0 : 3;
1816   Standard_Integer OrderU,OrderV;
1817   TColgp_Array2OfVec DerNUV(0,MaxOrder+3,0,MaxOrder+3);
1818   TColgp_Array2OfVec DerSurf(0,MaxOrder+4,0,MaxOrder+4);
1819   Standard_Real Umin,Umax,Vmin,Vmax;
1820   Bounds(Umin,Umax,Vmin,Vmax);
1821
1822   DerSurf.SetValue(1, 0, D1U);
1823   DerSurf.SetValue(0, 1, D1V);
1824   DerSurf.SetValue(1, 1, D2UV);
1825   DerSurf.SetValue(2, 0, D2U);
1826   DerSurf.SetValue(0, 2, D2V);
1827   DerSurf.SetValue(3, 0, D3U);
1828   DerSurf.SetValue(2, 1, D3UUV);
1829   DerSurf.SetValue(1, 2, D3UVV);
1830   DerSurf.SetValue(0, 3, D3V);
1831
1832
1833   //*********************
1834   Handle(Geom_BSplineSurface) L;
1835   Standard_Boolean IsOpposite=Standard_False;
1836   const Standard_Boolean AlongU = UOsculatingSurface(U,V,IsOpposite,L);
1837   const Standard_Boolean AlongV = VOsculatingSurface(U,V,IsOpposite,L);
1838   const Standard_Real signe = ((AlongV || AlongU) && IsOpposite)? -1. : 1.;
1839   derivatives(MaxOrder,3,U,V,basisSurf,3,3,AlongU,AlongV,L,DerNUV,DerSurf);
1840
1841   CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,NStatus,Normal,OrderU,OrderV);
1842   if (NStatus != CSLib_Defined) Geom_UndefinedValue::Raise();
1843
1844   P.SetXYZ(P.XYZ() + offsetValue * signe * Normal.XYZ());
1845
1846   D1U = DerSurf(1,0)
1847     + offsetValue * signe * CSLib::DNNormal(1,0,DerNUV,OrderU,OrderV);
1848   D1V = DerSurf(0,1)
1849     + offsetValue * signe * CSLib::DNNormal(0,1,DerNUV,OrderU,OrderV);
1850
1851   D2U = basisSurf->DN(U,V,2,0)
1852     + signe * offsetValue * CSLib::DNNormal(2,0,DerNUV,OrderU,OrderV);
1853   D2V = basisSurf->DN(U,V,0,2)
1854     + signe * offsetValue * CSLib::DNNormal(0,2,DerNUV,OrderU,OrderV);
1855   D2UV = basisSurf->DN(U,V,1,1)
1856     + signe * offsetValue * CSLib::DNNormal(1,1,DerNUV,OrderU,OrderV);
1857   D3U = basisSurf->DN(U,V,3,0)
1858     + signe * offsetValue * CSLib::DNNormal(3,0,DerNUV,OrderU,OrderV);
1859   D3V = basisSurf->DN(U,V,0,3)
1860     + signe * offsetValue * CSLib::DNNormal(0,3,DerNUV,OrderU,OrderV);
1861   D3UUV = basisSurf->DN(U,V,2,1)
1862     + signe * offsetValue * CSLib::DNNormal(2,1,DerNUV,OrderU,OrderV);
1863   D3UVV = basisSurf->DN(U,V,1,2)
1864     + signe * offsetValue * CSLib::DNNormal(1,2,DerNUV,OrderU,OrderV);
1865 }
1866
1867 //=======================================================================
1868 //function : SetDN
1869 //purpose  :
1870 //=======================================================================
1871
1872 gp_Vec Geom_OffsetSurface::SetDN ( const Standard_Real U, const Standard_Real V,
1873   const Standard_Integer Nu, const Standard_Integer Nv,
1874   const gp_Vec& D1U, const gp_Vec& D1V) const
1875 {
1876   const Standard_Real MagTol=0.000000001;
1877
1878   gp_Dir Normal;
1879   CSLib_NormalStatus NStatus;
1880   CSLib::Normal (D1U, D1V, MagTol, NStatus, Normal);
1881   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined)? 0 : 3;
1882   Standard_Integer OrderU,OrderV;
1883   TColgp_Array2OfVec DerNUV(0,MaxOrder+Nu,0,MaxOrder+Nu);
1884   TColgp_Array2OfVec DerSurf(0,MaxOrder+Nu+1,0,MaxOrder+Nv+1);
1885   Standard_Real Umin,Umax,Vmin,Vmax;
1886   Bounds(Umin,Umax,Vmin,Vmax);
1887
1888   DerSurf.SetValue(1, 0, D1U);
1889   DerSurf.SetValue(0, 1, D1V);
1890
1891   //*********************
1892   Handle(Geom_BSplineSurface) L;
1893   //   Is there any osculatingsurface along U or V;
1894   Standard_Boolean IsOpposite=Standard_False;
1895   const Standard_Boolean AlongU = UOsculatingSurface(U,V,IsOpposite,L);
1896   const Standard_Boolean AlongV = VOsculatingSurface(U,V,IsOpposite,L);
1897   const Standard_Real signe = ((AlongV || AlongU) && IsOpposite)? -1. : 1.;
1898   derivatives(MaxOrder,1,U,V,basisSurf,Nu,Nv,AlongU,AlongV,L,DerNUV,DerSurf);
1899
1900   CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,NStatus,Normal,OrderU,OrderV);
1901   if (NStatus != CSLib_Defined) Geom_UndefinedValue::Raise();
1902
1903   gp_Vec D = basisSurf->DN(U,V,Nu,Nv)
1904     + signe * offsetValue * CSLib::DNNormal(Nu,Nv,DerNUV,OrderU,OrderV);
1905
1906   return D;
1907 }