src/GeomEvaluator/GeomEvaluator_OffsetSurface.cxx

   1 // Created on: 2015-09-21
   2 // Copyright (c) 2015 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   6 // This library is free software; you can redistribute it and/or modify it under
   7 // the terms of the GNU Lesser General Public License version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 #include <GeomEvaluator_OffsetSurface.hxx>
  16
  17 #include <GeomAdaptor_HSurface.hxx>
  18 #include <CSLib.hxx>
  19 #include <CSLib_NormalStatus.hxx>
  20 #include <Geom_BezierSurface.hxx>
  21 #include <Geom_BSplineSurface.hxx>
  22 #include <Geom_UndefinedValue.hxx>
  23 #include <GeomAdaptor_HSurface.hxx>
  24 #include <gp_Vec.hxx>
  25 #include <Standard_RangeError.hxx>
  26
  27
  28 template<class SurfOrAdapt>
  29 static void derivatives(Standard_Integer theMaxOrder,
  30                         Standard_Integer theMinOrder,
  31                         const Standard_Real theU,
  32                         const Standard_Real theV,
  33                         const SurfOrAdapt& theBasisSurf,
  34                         const Standard_Integer theNU,
  35                         const Standard_Integer theNV,
  36                         const Standard_Boolean theAlongU,
  37                         const Standard_Boolean theAlongV,
  38                         const Handle(Geom_BSplineSurface)& theL,
  39                         TColgp_Array2OfVec& theDerNUV,
  40                         TColgp_Array2OfVec& theDerSurf);
  41
  42
  43
  44 GeomEvaluator_OffsetSurface::GeomEvaluator_OffsetSurface(
  45         const Handle(Geom_Surface)& theBase,
  46         const Standard_Real theOffset,
  47         const Handle(Geom_OsculatingSurface)& theOscSurf)
  48   : GeomEvaluator_Surface(),
  49     myBaseSurf(theBase),
  50     myOffset(theOffset),
  51     myOscSurf(theOscSurf)
  52 {
  53   if (!myOscSurf.IsNull())
  54     return; // osculating surface already exists
  55
  56   // Create osculating surface for B-spline and Besier surfaces only
  57   if (myBaseSurf->IsKind(STANDARD_TYPE(Geom_BSplineSurface)) ||
  58       myBaseSurf->IsKind(STANDARD_TYPE(Geom_BezierSurface)))
  59     myOscSurf = new Geom_OsculatingSurface(myBaseSurf, Precision::Confusion());
  60 }
  61
  62 GeomEvaluator_OffsetSurface::GeomEvaluator_OffsetSurface(
  63         const Handle(GeomAdaptor_HSurface)& theBase,
  64         const Standard_Real theOffset,
  65         const Handle(Geom_OsculatingSurface)& theOscSurf)
  66   : GeomEvaluator_Surface(),
  67     myBaseAdaptor(theBase),
  68     myOffset(theOffset),
  69     myOscSurf(theOscSurf)
  70 {
  71 }
  72
  73 void GeomEvaluator_OffsetSurface::D0(
  74     const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue) const
  75 {
  76   gp_Vec aD1U, aD1V;
  77   BaseD1(theU, theV, theValue, aD1U, aD1V);
  78   CalculateD0(theU, theV, theValue, aD1U, aD1V);
  79 }
  80
  81 void GeomEvaluator_OffsetSurface::D1(
  82     const Standard_Real theU, const Standard_Real theV,
  83     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V) const
  84 {
  85   gp_Vec aD2U, aD2V, aD2UV;
  86   BaseD2(theU, theV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
  87   CalculateD1(theU, theV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
  88 }
  89
  90 void GeomEvaluator_OffsetSurface::D2(
  91     const Standard_Real theU, const Standard_Real theV,
  92     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
  93     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV) const
  94 {
  95   gp_Vec aD3U, aD3V, aD3UUV, aD3UVV;
  96   BaseD3(theU, theV, theValue, theD1U, theD1V,
  97          theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
  98   CalculateD2(theU, theV, theValue, theD1U, theD1V,
  99               theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
 100 }
 101
 102 void GeomEvaluator_OffsetSurface::D3(
 103     const Standard_Real theU, const Standard_Real theV,
 104     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
 105     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
 106     gp_Vec& theD3U, gp_Vec& theD3V, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const
 107 {
 108   BaseD3(theU, theV, theValue, theD1U, theD1V,
 109          theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
 110   CalculateD3(theU, theV, theValue, theD1U, theD1V,
 111               theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
 112 }
 113
 114 gp_Vec GeomEvaluator_OffsetSurface::DN(
 115     const Standard_Real theU, const Standard_Real theV,
 116     const Standard_Integer theDerU, const Standard_Integer theDerV) const
 117 {
 118   Standard_RangeError_Raise_if(theDerU < 0, "GeomEvaluator_OffsetSurface::DN(): theDerU < 0");
 119   Standard_RangeError_Raise_if(theDerV < 0, "GeomEvaluator_OffsetSurface::DN(): theDerV < 0");
 120   Standard_RangeError_Raise_if(theDerU + theDerV < 1,
 121       "GeomEvaluator_OffsetSurface::DN(): theDerU + theDerV < 1");
 122
 123   gp_Pnt aP;
 124   gp_Vec aD1U, aD1V;
 125   BaseD1(theU, theV, aP, aD1U, aD1V);
 126   return CalculateDN(theU, theV, theDerU, theDerV, aD1U, aD1V);
 127 }
 128
 129
 130 void GeomEvaluator_OffsetSurface::BaseD0(const Standard_Real theU, const Standard_Real theV,
 131                                          gp_Pnt& theValue) const
 132 {
 133   if (!myBaseAdaptor.IsNull())
 134     myBaseAdaptor->D0(theU, theV, theValue);
 135   else
 136     myBaseSurf->D0(theU, theV, theValue);
 137 }
 138
 139 void GeomEvaluator_OffsetSurface::BaseD1(const Standard_Real theU, const Standard_Real theV,
 140                                          gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V) const
 141 {
 142   if (!myBaseAdaptor.IsNull())
 143     myBaseAdaptor->D1(theU, theV, theValue, theD1U, theD1V);
 144   else
 145     myBaseSurf->D1(theU, theV, theValue, theD1U, theD1V);
 146 }
 147
 148 void GeomEvaluator_OffsetSurface::BaseD2(const Standard_Real theU, const Standard_Real theV,
 149                                          gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
 150                                          gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV) const
 151 {
 152   if (!myBaseAdaptor.IsNull())
 153     myBaseAdaptor->D2(theU, theV, theValue, theD1U, theD1V, theD2U, theD2V, theD2UV);
 154   else
 155     myBaseSurf->D2(theU, theV, theValue, theD1U, theD1V, theD2U, theD2V, theD2UV);
 156 }
 157
 158 void GeomEvaluator_OffsetSurface::BaseD3(
 159     const Standard_Real theU, const Standard_Real theV,
 160     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
 161     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
 162     gp_Vec& theD3U, gp_Vec& theD3V, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const
 163 {
 164   if (!myBaseAdaptor.IsNull())
 165     myBaseAdaptor->D3(theU, theV, theValue, theD1U, theD1V,
 166     theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
 167   else
 168     myBaseSurf->D3(theU, theV, theValue, theD1U, theD1V,
 169     theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
 170 }
 171
 172
 173 void GeomEvaluator_OffsetSurface::CalculateD0(
 174     const Standard_Real theU, const Standard_Real theV,
 175     gp_Pnt& theValue,
 176     const gp_Vec& theD1U, const gp_Vec& theD1V) const
 177 {
 178   const Standard_Real MagTol = 1.e-9;
 179
 180   // Normalize derivatives before normal calculation because it gives more stable result.
 181   // There will be normalized only derivatives greater than 1.0 to avoid differences in last significant digit
 182   gp_Vec aD1U(theD1U);
 183   gp_Vec aD1V(theD1V);
 184   Standard_Real aD1UNorm2 = aD1U.SquareMagnitude();
 185   Standard_Real aD1VNorm2 = aD1V.SquareMagnitude();
 186   if (aD1UNorm2 > 1.0)
 187     aD1U /= Sqrt(aD1UNorm2);
 188   if (aD1VNorm2 > 1.0)
 189     aD1V /= Sqrt(aD1VNorm2);
 190
 191   gp_Vec aNorm = aD1U.Crossed(aD1V);
 192   if (aNorm.SquareMagnitude() > MagTol * MagTol)
 193   {
 194     // Non singular case. Simple computations.
 195     aNorm.Normalize();
 196     theValue.SetXYZ(theValue.XYZ() + myOffset * aNorm.XYZ());
 197   }
 198   else
 199   {
 200     const Standard_Integer MaxOrder = 3;
 201
 202     Handle(Geom_BSplineSurface) L;
 203     Standard_Boolean isOpposite = Standard_False;
 204     Standard_Boolean AlongU = Standard_False;
 205     Standard_Boolean AlongV = Standard_False;
 206     if (!myOscSurf.IsNull())
 207     {
 208       AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
 209       AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
 210     }
 211     const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
 212
 213     TColgp_Array2OfVec DerNUV(0, MaxOrder, 0, MaxOrder);
 214     TColgp_Array2OfVec DerSurf(0, MaxOrder + 1, 0, MaxOrder + 1);
 215     Standard_Integer OrderU, OrderV;
 216     Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
 217     Bounds(Umin, Umax, Vmin, Vmax);
 218
 219     DerSurf.SetValue(1, 0, theD1U);
 220     DerSurf.SetValue(0, 1, theD1V);
 221     if (!myBaseSurf.IsNull())
 222       derivatives(MaxOrder, 1, theU, theV, myBaseSurf, 0, 0, AlongU, AlongV, L, DerNUV, DerSurf);
 223     else
 224       derivatives(MaxOrder, 1, theU, theV, myBaseAdaptor, 0, 0, AlongU, AlongV, L, DerNUV, DerSurf);
 225
 226     gp_Dir Normal;
 227     CSLib_NormalStatus NStatus;
 228     CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
 229                   NStatus, Normal, OrderU, OrderV);
 230     if (NStatus == CSLib_InfinityOfSolutions)
 231     {
 232       // Replace zero derivative and try to calculate normal
 233       gp_Vec aNewDU = theD1U;
 234       gp_Vec aNewDV = theD1V;
 235       if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, MagTol * MagTol))
 236         CSLib::Normal(aNewDU, aNewDV, MagTol, NStatus, Normal);
 237     }
 238
 239     if (NStatus != CSLib_Defined)
 240       Geom_UndefinedValue::Raise(
 241           "GeomEvaluator_OffsetSurface::CalculateD0(): Unable to calculate normal");
 242
 243     theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
 244   }
 245 }
 246
 247 void GeomEvaluator_OffsetSurface::CalculateD1(
 248     const Standard_Real theU, const Standard_Real theV,
 249     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
 250     const gp_Vec& theD2U, const gp_Vec& theD2V, const gp_Vec& theD2UV) const
 251 {
 252   const Standard_Real MagTol = 1.e-9;
 253
 254   // Check offset side.
 255   Handle(Geom_BSplineSurface) L;
 256   Standard_Boolean isOpposite = Standard_False;
 257   Standard_Boolean AlongU = Standard_False;
 258   Standard_Boolean AlongV = Standard_False;
 259   if (!myOscSurf.IsNull())
 260   {
 261     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
 262     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
 263   }
 264   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
 265
 266   // Normalize derivatives before normal calculation because it gives more stable result.
 267   // There will be normalized only derivatives greater than 1.0 to avoid differences in last significant digit
 268   gp_Vec aD1U(theD1U);
 269   gp_Vec aD1V(theD1V);
 270   Standard_Real aD1UNorm2 = aD1U.SquareMagnitude();
 271   Standard_Real aD1VNorm2 = aD1V.SquareMagnitude();
 272   if (aD1UNorm2 > 1.0)
 273     aD1U /= Sqrt(aD1UNorm2);
 274   if (aD1VNorm2 > 1.0)
 275     aD1V /= Sqrt(aD1VNorm2);
 276
 277   Standard_Integer MaxOrder = 0;
 278   gp_Vec aNorm = aD1U.Crossed(aD1V);
 279   if (aNorm.SquareMagnitude() > MagTol * MagTol)
 280   {
 281     if (!AlongV && !AlongU)
 282     {
 283       // AlongU or AlongV leads to more complex D1 computation
 284       // Try to compute D0 and D1 much simpler
 285       aNorm.Normalize();
 286       theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * aNorm.XYZ());
 287
 288       gp_Vec aN0(aNorm.XYZ()), aN1U, aN1V;
 289       Standard_Real aScale = (theD1U^theD1V).Dot(aN0);
 290       aN1U.SetX(theD2U.Y() * theD1V.Z() + theD1U.Y() * theD2UV.Z()
 291         - theD2U.Z() * theD1V.Y() - theD1U.Z() * theD2UV.Y());
 292       aN1U.SetY((theD2U.X() * theD1V.Z() + theD1U.X() * theD2UV.Z()
 293         - theD2U.Z() * theD1V.X() - theD1U.Z() * theD2UV.X()) * -1.0);
 294       aN1U.SetZ(theD2U.X() * theD1V.Y() + theD1U.X() * theD2UV.Y()
 295         - theD2U.Y() * theD1V.X() - theD1U.Y() * theD2UV.X());
 296       Standard_Real aScaleU = aN1U.Dot(aN0);
 297       aN1U.Subtract(aScaleU * aN0);
 298       aN1U /= aScale;
 299
 300       aN1V.SetX(theD2UV.Y() * theD1V.Z() + theD2V.Z() * theD1U.Y()
 301         - theD2UV.Z() * theD1V.Y() - theD2V.Y() * theD1U.Z());
 302       aN1V.SetY((theD2UV.X() * theD1V.Z() + theD2V.Z() * theD1U.X()
 303         - theD2UV.Z() * theD1V.X() - theD2V.X() * theD1U.Z()) * -1.0);
 304       aN1V.SetZ(theD2UV.X() * theD1V.Y() + theD2V.Y() * theD1U.X()
 305         - theD2UV.Y() * theD1V.X() - theD2V.X() * theD1U.Y());
 306       Standard_Real aScaleV = aN1V.Dot(aN0);
 307       aN1V.Subtract(aScaleV * aN0);
 308       aN1V /= aScale;
 309
 310       theD1U += myOffset * aSign * aN1U;
 311       theD1V += myOffset * aSign * aN1V;
 312
 313       return;
 314     }
 315   }
 316   else
 317     MaxOrder = 3;
 318
 319   Standard_Integer OrderU, OrderV;
 320   TColgp_Array2OfVec DerNUV(0, MaxOrder + 1, 0, MaxOrder + 1);
 321   TColgp_Array2OfVec DerSurf(0, MaxOrder + 2, 0, MaxOrder + 2);
 322   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
 323   Bounds(Umin, Umax, Vmin, Vmax);
 324
 325   DerSurf.SetValue(1, 0, theD1U);
 326   DerSurf.SetValue(0, 1, theD1V);
 327   DerSurf.SetValue(1, 1, theD2UV);
 328   DerSurf.SetValue(2, 0, theD2U);
 329   DerSurf.SetValue(0, 2, theD2V);
 330   if (!myBaseSurf.IsNull())
 331     derivatives(MaxOrder, 2, theU, theV, myBaseSurf, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
 332   else
 333     derivatives(MaxOrder, 2, theU, theV, myBaseAdaptor, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
 334
 335   gp_Dir Normal;
 336   CSLib_NormalStatus NStatus;
 337   CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax, NStatus, Normal, OrderU, OrderV);
 338   if (NStatus == CSLib_InfinityOfSolutions)
 339   {
 340     gp_Vec aNewDU = theD1U;
 341     gp_Vec aNewDV = theD1V;
 342     // Replace zero derivative and try to calculate normal
 343     if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, MagTol * MagTol))
 344     {
 345       DerSurf.SetValue(1, 0, aNewDU);
 346       DerSurf.SetValue(0, 1, aNewDV);
 347       if (!myBaseSurf.IsNull())
 348         derivatives(MaxOrder, 2, theU, theV, myBaseSurf, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
 349       else
 350         derivatives(MaxOrder, 2, theU, theV, myBaseAdaptor, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
 351       CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
 352                     NStatus, Normal, OrderU, OrderV);
 353     }
 354   }
 355
 356   if (NStatus != CSLib_Defined)
 357     Geom_UndefinedValue::Raise(
 358         "GeomEvaluator_OffsetSurface::CalculateD1(): Unable to calculate normal");
 359
 360   theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
 361
 362   theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV);
 363   theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV);
 364 }
 365
 366 void GeomEvaluator_OffsetSurface::CalculateD2(
 367     const Standard_Real theU, const Standard_Real theV,
 368     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
 369     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
 370     const gp_Vec& theD3U, const gp_Vec& theD3V, const gp_Vec& theD3UUV, const gp_Vec& theD3UVV) const
 371 {
 372   const Standard_Real MagTol = 1.e-9;
 373
 374   gp_Dir Normal;
 375   CSLib_NormalStatus NStatus;
 376   CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
 377
 378   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
 379   Standard_Integer OrderU, OrderV;
 380   TColgp_Array2OfVec DerNUV(0, MaxOrder + 2, 0, MaxOrder + 2);
 381   TColgp_Array2OfVec DerSurf(0, MaxOrder + 3, 0, MaxOrder + 3);
 382
 383   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
 384   Bounds(Umin, Umax, Vmin, Vmax);
 385
 386   DerSurf.SetValue(1, 0, theD1U);
 387   DerSurf.SetValue(0, 1, theD1V);
 388   DerSurf.SetValue(1, 1, theD2UV);
 389   DerSurf.SetValue(2, 0, theD2U);
 390   DerSurf.SetValue(0, 2, theD2V);
 391   DerSurf.SetValue(3, 0, theD3U);
 392   DerSurf.SetValue(2, 1, theD3UUV);
 393   DerSurf.SetValue(1, 2, theD3UVV);
 394   DerSurf.SetValue(0, 3, theD3V);
 395   //*********************
 396
 397   Handle(Geom_BSplineSurface) L;
 398   Standard_Boolean isOpposite = Standard_False;
 399   Standard_Boolean AlongU = Standard_False;
 400   Standard_Boolean AlongV = Standard_False;
 401   if (!myOscSurf.IsNull())
 402   {
 403     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
 404     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
 405   }
 406   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
 407
 408   if (!myBaseSurf.IsNull())
 409     derivatives(MaxOrder, 3, theU, theV, myBaseSurf, 2, 2, AlongU, AlongV, L, DerNUV, DerSurf);
 410   else
 411     derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 2, 2, AlongU, AlongV, L, DerNUV, DerSurf);
 412
 413   CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
 414                 NStatus, Normal, OrderU, OrderV);
 415   if (NStatus != CSLib_Defined)
 416     Geom_UndefinedValue::Raise(
 417         "GeomEvaluator_OffsetSurface::CalculateD2(): Unable to calculate normal");
 418
 419   theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
 420
 421   theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV);
 422   theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV);
 423
 424   if (!myBaseSurf.IsNull())
 425   {
 426     theD2U = myBaseSurf->DN(theU, theV, 2, 0);
 427     theD2V = myBaseSurf->DN(theU, theV, 0, 2);
 428     theD2UV = myBaseSurf->DN(theU, theV, 1, 1);
 429   }
 430   else
 431   {
 432     theD2U = myBaseAdaptor->DN(theU, theV, 2, 0);
 433     theD2V = myBaseAdaptor->DN(theU, theV, 0, 2);
 434     theD2UV = myBaseAdaptor->DN(theU, theV, 1, 1);
 435   }
 436
 437   theD2U += aSign * myOffset * CSLib::DNNormal(2, 0, DerNUV, OrderU, OrderV);
 438   theD2V += aSign * myOffset * CSLib::DNNormal(0, 2, DerNUV, OrderU, OrderV);
 439   theD2UV += aSign * myOffset * CSLib::DNNormal(1, 1, DerNUV, OrderU, OrderV);
 440 }
 441
 442 void GeomEvaluator_OffsetSurface::CalculateD3(
 443     const Standard_Real theU, const Standard_Real theV,
 444     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
 445     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
 446     gp_Vec& theD3U, gp_Vec& theD3V, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const
 447 {
 448   const Standard_Real MagTol = 1.e-9;
 449
 450   gp_Dir Normal;
 451   CSLib_NormalStatus NStatus;
 452   CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
 453   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
 454   Standard_Integer OrderU, OrderV;
 455   TColgp_Array2OfVec DerNUV(0, MaxOrder + 3, 0, MaxOrder + 3);
 456   TColgp_Array2OfVec DerSurf(0, MaxOrder + 4, 0, MaxOrder + 4);
 457   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
 458   Bounds(Umin, Umax, Vmin, Vmax);
 459
 460   DerSurf.SetValue(1, 0, theD1U);
 461   DerSurf.SetValue(0, 1, theD1V);
 462   DerSurf.SetValue(1, 1, theD2UV);
 463   DerSurf.SetValue(2, 0, theD2U);
 464   DerSurf.SetValue(0, 2, theD2V);
 465   DerSurf.SetValue(3, 0, theD3U);
 466   DerSurf.SetValue(2, 1, theD3UUV);
 467   DerSurf.SetValue(1, 2, theD3UVV);
 468   DerSurf.SetValue(0, 3, theD3V);
 469
 470
 471   //*********************
 472   Handle(Geom_BSplineSurface) L;
 473   Standard_Boolean isOpposite = Standard_False;
 474   Standard_Boolean AlongU = Standard_False;
 475   Standard_Boolean AlongV = Standard_False;
 476   if (!myOscSurf.IsNull())
 477   {
 478     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
 479     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
 480   }
 481   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
 482
 483   if (!myBaseSurf.IsNull())
 484     derivatives(MaxOrder, 3, theU, theV, myBaseSurf, 3, 3, AlongU, AlongV, L, DerNUV, DerSurf);
 485   else
 486     derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 3, 3, AlongU, AlongV, L, DerNUV, DerSurf);
 487
 488   CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
 489                 NStatus, Normal, OrderU, OrderV);
 490   if (NStatus != CSLib_Defined)
 491     Geom_UndefinedValue::Raise(
 492         "GeomEvaluator_OffsetSurface::CalculateD3(): Unable to calculate normal");
 493
 494   theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
 495
 496   theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV);
 497   theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV);
 498
 499   if (!myBaseSurf.IsNull())
 500   {
 501     theD2U   = myBaseSurf->DN(theU, theV, 2, 0);
 502     theD2V   = myBaseSurf->DN(theU, theV, 0, 2);
 503     theD2UV  = myBaseSurf->DN(theU, theV, 1, 1);
 504     theD3U   = myBaseSurf->DN(theU, theV, 3, 0);
 505     theD3V   = myBaseSurf->DN(theU, theV, 0, 3);
 506     theD3UUV = myBaseSurf->DN(theU, theV, 2, 1);
 507     theD3UVV = myBaseSurf->DN(theU, theV, 1, 2);
 508   }
 509   else
 510   {
 511     theD2U   = myBaseAdaptor->DN(theU, theV, 2, 0);
 512     theD2V   = myBaseAdaptor->DN(theU, theV, 0, 2);
 513     theD2UV  = myBaseAdaptor->DN(theU, theV, 1, 1);
 514     theD3U   = myBaseAdaptor->DN(theU, theV, 3, 0);
 515     theD3V   = myBaseAdaptor->DN(theU, theV, 0, 3);
 516     theD3UUV = myBaseAdaptor->DN(theU, theV, 2, 1);
 517     theD3UVV = myBaseAdaptor->DN(theU, theV, 1, 2);
 518   }
 519
 520   theD2U   += aSign * myOffset * CSLib::DNNormal(2, 0, DerNUV, OrderU, OrderV);
 521   theD2V   += aSign * myOffset * CSLib::DNNormal(0, 2, DerNUV, OrderU, OrderV);
 522   theD2UV  += aSign * myOffset * CSLib::DNNormal(1, 1, DerNUV, OrderU, OrderV);
 523   theD3U   += aSign * myOffset * CSLib::DNNormal(3, 0, DerNUV, OrderU, OrderV);
 524   theD3V   += aSign * myOffset * CSLib::DNNormal(0, 3, DerNUV, OrderU, OrderV);
 525   theD3UUV += aSign * myOffset * CSLib::DNNormal(2, 1, DerNUV, OrderU, OrderV);
 526   theD3UVV += aSign * myOffset * CSLib::DNNormal(1, 2, DerNUV, OrderU, OrderV);
 527 }
 528
 529 gp_Vec GeomEvaluator_OffsetSurface::CalculateDN(
 530     const Standard_Real theU, const Standard_Real theV,
 531     const Standard_Integer theNu, const Standard_Integer theNv,
 532     const gp_Vec& theD1U, const gp_Vec& theD1V) const
 533 {
 534   const Standard_Real MagTol = 1.e-9;
 535
 536   gp_Dir Normal;
 537   CSLib_NormalStatus NStatus;
 538   CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
 539   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
 540   Standard_Integer OrderU, OrderV;
 541   TColgp_Array2OfVec DerNUV(0, MaxOrder + theNu, 0, MaxOrder + theNu);
 542   TColgp_Array2OfVec DerSurf(0, MaxOrder + theNu + 1, 0, MaxOrder + theNv + 1);
 543
 544   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
 545   Bounds(Umin, Umax, Vmin, Vmax);
 546
 547   DerSurf.SetValue(1, 0, theD1U);
 548   DerSurf.SetValue(0, 1, theD1V);
 549
 550   //*********************
 551   Handle(Geom_BSplineSurface) L;
 552   //   Is there any osculatingsurface along U or V;
 553   Standard_Boolean isOpposite = Standard_False;
 554   Standard_Boolean AlongU = Standard_False;
 555   Standard_Boolean AlongV = Standard_False;
 556   if (!myOscSurf.IsNull())
 557   {
 558     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
 559     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
 560   }
 561   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
 562
 563   if (!myBaseSurf.IsNull())
 564     derivatives(MaxOrder, 1, theU, theV, myBaseSurf, theNu, theNv, AlongU, AlongV, L, DerNUV, DerSurf);
 565   else
 566     derivatives(MaxOrder, 1, theU, theV, myBaseAdaptor, theNu, theNv, AlongU, AlongV, L, DerNUV, DerSurf);
 567
 568   CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
 569                 NStatus, Normal, OrderU, OrderV);
 570   if (NStatus != CSLib_Defined)
 571     Geom_UndefinedValue::Raise(
 572         "GeomEvaluator_OffsetSurface::CalculateDN(): Unable to calculate normal");
 573
 574   gp_Vec D;
 575   if (!myBaseSurf.IsNull())
 576     D = myBaseSurf->DN(theU, theV, theNu, theNv);
 577   else
 578     D = myBaseAdaptor->DN(theU, theV, theNu, theNv);
 579
 580   D += aSign * myOffset * CSLib::DNNormal(theNu, theNv, DerNUV, OrderU, OrderV);
 581   return D;
 582 }
 583
 584
 585 void GeomEvaluator_OffsetSurface::Bounds(Standard_Real& theUMin, Standard_Real& theUMax,
 586                                          Standard_Real& theVMin, Standard_Real& theVMax) const
 587 {
 588   if (!myBaseSurf.IsNull())
 589     myBaseSurf->Bounds(theUMin, theUMax, theVMin, theVMax);
 590   else
 591   {
 592     theUMin = myBaseAdaptor->FirstUParameter();
 593     theUMax = myBaseAdaptor->LastUParameter();
 594     theVMin = myBaseAdaptor->FirstVParameter();
 595     theVMax = myBaseAdaptor->LastVParameter();
 596   }
 597 }
 598
 599
 600 Standard_Boolean GeomEvaluator_OffsetSurface::ReplaceDerivative(
 601     const Standard_Real theU, const Standard_Real theV,
 602     gp_Vec& theDU, gp_Vec& theDV,
 603     const Standard_Real theSquareTol) const
 604 {
 605   Standard_Boolean isReplaceDU = theDU.SquareMagnitude() < theSquareTol;
 606   Standard_Boolean isReplaceDV = theDV.SquareMagnitude() < theSquareTol;
 607   Standard_Boolean isReplaced = Standard_False;
 608   if (isReplaceDU ^ isReplaceDV)
 609   {
 610     Standard_Real aU = theU;
 611     Standard_Real aV = theV;
 612     Standard_Real aUMin = 0, aUMax = 0, aVMin = 0, aVMax = 0;
 613     Bounds(aUMin, aUMax, aVMin, aVMax);
 614
 615     // Calculate step along non-zero derivative
 616     Standard_Real aStep;
 617     Handle(Adaptor3d_HSurface) aSurfAdapt;
 618     if (!myBaseAdaptor.IsNull())
 619       aSurfAdapt = myBaseAdaptor;
 620     else
 621       aSurfAdapt = new GeomAdaptor_HSurface(myBaseSurf);
 622     if (isReplaceDV)
 623     {
 624       aStep = Precision::Confusion() * theDU.Magnitude();
 625       if (aStep > aUMax - aUMin)
 626         aStep = (aUMax - aUMin) / 100.;
 627     }
 628     else
 629     {
 630       aStep = Precision::Confusion() * theDV.Magnitude();
 631       if (aStep > aVMax - aVMin)
 632         aStep = (aVMax - aVMin) / 100.;
 633     }
 634
 635     gp_Pnt aP;
 636     gp_Vec aDU, aDV;
 637     // Step away from currect parametric coordinates and calculate derivatives once again.
 638     // Replace zero derivative by the obtained.
 639     for (Standard_Real aStepSign = -1.0; aStepSign <= 1.0 && !isReplaced; aStepSign += 2.0)
 640     {
 641       if (isReplaceDV)
 642       {
 643         aU = theU + aStepSign * aStep;
 644         if (aU < aUMin || aU > aUMax)
 645           continue;
 646       }
 647       else
 648       {
 649         aV = theV + aStepSign * aStep;
 650         if (aV < aVMin || aV > aVMax)
 651           continue;
 652       }
 653
 654       BaseD1(aU, aV, aP, aDU, aDV);
 655
 656       if (isReplaceDU && aDU.SquareMagnitude() > theSquareTol)
 657       {
 658         theDU = aDU;
 659         isReplaced = Standard_True;
 660       }
 661       if (isReplaceDV && aDV.SquareMagnitude() > theSquareTol)
 662       {
 663         theDV = aDV;
 664         isReplaced = Standard_True;
 665       }
 666     }
 667   }
 668   return isReplaced;
 669 }
 670
 671
 672
 673
 674
 675 // =====================   Auxiliary functions   ===================================
 676 template<class SurfOrAdapt>
 677 void derivatives(Standard_Integer theMaxOrder,
 678                  Standard_Integer theMinOrder,
 679                  const Standard_Real theU,
 680                  const Standard_Real theV,
 681                  const SurfOrAdapt& theBasisSurf,
 682                  const Standard_Integer theNU,
 683                  const Standard_Integer theNV,
 684                  const Standard_Boolean theAlongU,
 685                  const Standard_Boolean theAlongV,
 686                  const Handle(Geom_BSplineSurface)& theL,
 687                  TColgp_Array2OfVec& theDerNUV,
 688                  TColgp_Array2OfVec& theDerSurf)
 689 {
 690   Standard_Integer i, j;
 691   gp_Pnt P;
 692   gp_Vec DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3UUV, DL3UVV, DL3V;
 693
 694   if (theAlongU || theAlongV)
 695   {
 696     theMaxOrder = 0;
 697     TColgp_Array2OfVec DerSurfL(0, theMaxOrder + theNU + 1, 0, theMaxOrder + theNV + 1);
 698     switch (theMinOrder)
 699     {
 700     case 1:
 701       theL->D1(theU, theV, P, DL1U, DL1V);
 702       DerSurfL.SetValue(1, 0, DL1U);
 703       DerSurfL.SetValue(0, 1, DL1V);
 704       break;
 705     case 2:
 706       theL->D2(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV);
 707       DerSurfL.SetValue(1, 0, DL1U);
 708       DerSurfL.SetValue(0, 1, DL1V);
 709       DerSurfL.SetValue(1, 1, DL2UV);
 710       DerSurfL.SetValue(2, 0, DL2U);
 711       DerSurfL.SetValue(0, 2, DL2V);
 712       break;
 713     case 3:
 714       theL->D3(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3V, DL3UUV, DL3UVV);
 715       DerSurfL.SetValue(1, 0, DL1U);
 716       DerSurfL.SetValue(0, 1, DL1V);
 717       DerSurfL.SetValue(1, 1, DL2UV);
 718       DerSurfL.SetValue(2, 0, DL2U);
 719       DerSurfL.SetValue(0, 2, DL2V);
 720       DerSurfL.SetValue(3, 0, DL3U);
 721       DerSurfL.SetValue(2, 1, DL3UUV);
 722       DerSurfL.SetValue(1, 2, DL3UVV);
 723       DerSurfL.SetValue(0, 3, DL3V);
 724       break;
 725     default:
 726       break;
 727     }
 728
 729     if (theNU <= theNV)
 730     {
 731       for (i = 0; i <= theMaxOrder + 1 + theNU; i++)
 732         for (j = i; j <= theMaxOrder + theNV + 1; j++)
 733           if (i + j > theMinOrder)
 734           {
 735             DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
 736             theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
 737             if (i != j && j <= theNU + 1)
 738             {
 739               theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
 740               DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
 741             }
 742           }
 743     }
 744     else
 745     {
 746       for (j = 0; j <= theMaxOrder + 1 + theNV; j++)
 747         for (i = j; i <= theMaxOrder + theNU + 1; i++)
 748           if (i + j > theMinOrder)
 749           {
 750             DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
 751             theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
 752             if (i != j && i <= theNV + 1)
 753             {
 754               theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
 755               DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
 756             }
 757           }
 758     }
 759     for (i = 0; i <= theMaxOrder + theNU; i++)
 760       for (j = 0; j <= theMaxOrder + theNV; j++)
 761       {
 762         if (theAlongU)
 763           theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, DerSurfL, theDerSurf));
 764         if (theAlongV)
 765           theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf, DerSurfL));
 766       }
 767   }
 768   else
 769   {
 770     for (i = 0; i <= theMaxOrder + theNU+ 1; i++)
 771       for (j = i; j <= theMaxOrder + theNV + 1; j++)
 772         if (i + j > theMinOrder)
 773         {
 774           theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
 775           if (i != j)
 776             theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
 777         }
 778     for (i = 0; i <= theMaxOrder + theNU; i++)
 779       for (j = 0; j <= theMaxOrder + theNV; j++)
 780         theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf));
 781   }
 782 }
 783