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