src/CSLib/CSLib_Offset.cxx

   1 // Copyright (c) 2015-... OPEN CASCADE SAS
   2 //
   3 // This file is part of Open CASCADE Technology software library.
   4 //
   5 // This library is free software; you can redistribute it and/or modify it under
   6 // the terms of the GNU Lesser General Public License version 2.1 as published
   7 // by the Free Software Foundation, with special exception defined in the file
   8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
   9 // distribution for complete text of the license and disclaimer of any warranty.
  10 //
  11 // Alternatively, this file may be used under the terms of Open CASCADE
  12 // commercial license or contractual agreement.
  13
  14 #include <CSLib_Offset.hxx>
  15
  16 #include <gp_Dir2d.hxx>
  17 #include <gp_XY.hxx>
  18
  19 #include <Standard_NullValue.hxx>
  20
  21
  22 // ==========  Offset values for 2D curves  ==========
  23
  24 void CSLib_Offset::D0(const gp_Pnt2d&  theBasePoint,
  25                       const gp_Vec2d&  theBaseDeriv,
  26                       Standard_Real    theOffset,
  27                       Standard_Boolean , // unused
  28                       gp_Pnt2d&        theResPoint)
  29 {
  30   if (theBaseDeriv.SquareMagnitude() <= gp::Resolution())
  31     Standard_NullValue::Raise("CSLib_Offset: Undefined normal vector "
  32                               "because tangent vector has zero-magnitude!");
  33
  34   gp_Dir2d aNormal(theBaseDeriv.Y(), -theBaseDeriv.X());
  35   theResPoint.SetCoord(theBasePoint.X() + aNormal.X() * theOffset,
  36                        theBasePoint.Y() + aNormal.Y() * theOffset);
  37 }
  38
  39 void CSLib_Offset::D1(const gp_Pnt2d& theBasePoint,
  40                       const gp_Vec2d& theBaseD1,
  41                       const gp_Vec2d& theBaseD2,
  42                       Standard_Real   theOffset,
  43                       Standard_Boolean , // unused
  44                       gp_Pnt2d&       theResPoint,
  45                       gp_Vec2d&       theResDeriv)
  46 {
  47   // P(u) = p(u) + Offset * Ndir / R
  48   // with R = || p' ^ Z|| and Ndir = P' ^ Z
  49
  50   // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
  51
  52   gp_XY Ndir(theBaseD1.Y(), -theBaseD1.X());
  53   gp_XY DNdir(theBaseD2.Y(), -theBaseD2.X());
  54   Standard_Real R2 = Ndir.SquareModulus();
  55   Standard_Real R  = Sqrt (R2);
  56   Standard_Real R3 = R * R2;
  57   Standard_Real Dr = Ndir.Dot(DNdir);
  58   if (R3 <= gp::Resolution())
  59   {
  60     if (R2 <= gp::Resolution())
  61       Standard_NullValue::Raise("CSLib_Offset: Null derivative");
  62     //We try another computation but the stability is not very good.
  63     DNdir.Multiply(R);
  64     DNdir.Subtract(Ndir.Multiplied(Dr / R));
  65     DNdir.Multiply(theOffset / R2);
  66   }
  67   else
  68   {
  69     // Same computation as IICURV in EUCLID-IS because the stability is better
  70     DNdir.Multiply(theOffset / R);
  71     DNdir.Subtract(Ndir.Multiplied(theOffset * Dr / R3));
  72   }
  73
  74   // P(u)
  75   D0(theBasePoint, theBaseD1, theOffset, Standard_False, theResPoint);
  76   // P'(u)
  77   theResDeriv = theBaseD1.Added(gp_Vec2d(DNdir));
  78 }
  79
  80 void CSLib_Offset::D2(const gp_Pnt2d&  theBasePoint,
  81                       const gp_Vec2d&  theBaseD1,
  82                       const gp_Vec2d&  theBaseD2,
  83                       const gp_Vec2d&  theBaseD3,
  84                       Standard_Real    theOffset,
  85                       Standard_Boolean theIsDirectionChange,
  86                       gp_Pnt2d&        theResPoint,
  87                       gp_Vec2d&        theResD1,
  88                       gp_Vec2d&        theResD2)
  89 {
  90   // P(u) = p(u) + Offset * Ndir / R
  91   // with R = || p' ^ Z|| and Ndir = P' ^ Z
  92
  93   // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
  94
  95   // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
  96   //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
  97
  98   gp_XY Ndir(theBaseD1.Y(), -theBaseD1.X());
  99   gp_XY DNdir(theBaseD2.Y(), -theBaseD2.X());
 100   gp_XY D2Ndir(theBaseD3.Y(), -theBaseD3.X());
 101   Standard_Real R2  = Ndir.SquareModulus();
 102   Standard_Real R   = Sqrt(R2);
 103   Standard_Real R3  = R2 * R;
 104   Standard_Real R4  = R2 * R2;
 105   Standard_Real R5  = R3 * R2;
 106   Standard_Real Dr  = Ndir.Dot(DNdir);
 107   Standard_Real D2r = Ndir.Dot(D2Ndir) + DNdir.Dot (DNdir);
 108   if (R5 <= gp::Resolution())
 109   {
 110     if (R4 <= gp::Resolution())
 111       Standard_NullValue::Raise("CSLib_Offset: Null derivative");
 112     //We try another computation but the stability is not very good dixit ISG.
 113     // V2 = P" (U) :
 114     D2Ndir.Subtract(DNdir.Multiplied (2.0 * Dr / R2));
 115     D2Ndir.Add(Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
 116     D2Ndir.Multiply(theOffset / R);
 117
 118     // V1 = P' (U) :
 119     DNdir.Multiply(R);
 120     DNdir.Subtract(Ndir.Multiplied(Dr / R));
 121     DNdir.Multiply(theOffset / R2);
 122   }
 123   else
 124   {
 125     // Same computation as IICURV in EUCLID-IS because the stability is better.
 126     // V2 = P" (U) :
 127     D2Ndir.Multiply(theOffset / R);
 128     D2Ndir.Subtract(DNdir.Multiplied (2.0 * theOffset * Dr / R3));
 129     D2Ndir.Add (Ndir.Multiplied(theOffset * (((3.0 * Dr * Dr) / R5) - (D2r / R3))));
 130
 131     // V1 = P' (U)
 132     DNdir.Multiply(theOffset / R);
 133     DNdir.Subtract(Ndir.Multiplied(theOffset * Dr / R3));
 134   }
 135
 136   // P(u) :
 137   D0(theBasePoint, theBaseD1, theOffset, theIsDirectionChange, theResPoint);
 138   // P'(u) :
 139   theResD1 = theBaseD1.Added(gp_Vec2d(DNdir));
 140   // P"(u) :
 141   if (theIsDirectionChange)
 142     theResD2 = -theBaseD2;
 143   else
 144     theResD2 = theBaseD2;
 145   theResD2.Add(gp_Vec2d(D2Ndir));
 146 }
 147
 148 void CSLib_Offset::D3(const gp_Pnt2d&  theBasePoint,
 149                       const gp_Vec2d&  theBaseD1,
 150                       const gp_Vec2d&  theBaseD2,
 151                       const gp_Vec2d&  theBaseD3,
 152                       const gp_Vec2d&  theBaseD4,
 153                       Standard_Real    theOffset,
 154                       Standard_Boolean theIsDirectionChange,
 155                       gp_Pnt2d&        theResPoint,
 156                       gp_Vec2d&        theResD1,
 157                       gp_Vec2d&        theResD2,
 158                       gp_Vec2d&        theResD3)
 159 {
 160   // P(u) = p(u) + Offset * Ndir / R
 161   // with R = || p' ^ Z|| and Ndir = P' ^ Z
 162
 163   // P'(u)  = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 164
 165   // P"(u)  = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
 166   //          Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 167
 168   // P"'(u) = p"'(u) + (Offset / R) * (D3Ndir - (3.0 * Dr/R**2 ) * D2Ndir -
 169   //          (3.0 * D2r / R2) * DNdir) + (3.0 * Dr * Dr / R4) * DNdir -
 170   //          (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir +
 171   //          (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir
 172
 173   gp_XY Ndir(theBaseD1.Y(), -theBaseD1.X());
 174   gp_XY DNdir(theBaseD2.Y(), -theBaseD2.X());
 175   gp_XY D2Ndir(theBaseD3.Y(), -theBaseD3.X());
 176   gp_XY D3Ndir(theBaseD4.Y(), -theBaseD4.X());
 177   Standard_Real R2  = Ndir.SquareModulus();
 178   Standard_Real R   = Sqrt (R2);
 179   Standard_Real R3  = R2 * R;
 180   Standard_Real R4  = R2 * R2;
 181   Standard_Real R5  = R3 * R2;
 182   Standard_Real R6  = R3 * R3;
 183   Standard_Real R7  = R5 * R2;
 184   Standard_Real Dr  = Ndir.Dot(DNdir);
 185   Standard_Real D2r = Ndir.Dot(D2Ndir) + DNdir.Dot (DNdir);
 186   Standard_Real D3r = Ndir.Dot(D3Ndir) + 3.0 * DNdir.Dot (D2Ndir);
 187
 188   if (R7 <= gp::Resolution())
 189   {
 190     if (R6 <= gp::Resolution())
 191       Standard_NullValue::Raise("CSLib_Offset: Null derivative");
 192     //We try another computation but the stability is not very good dixit ISG.
 193     // V3 = P"' (U) :
 194     D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * theOffset * Dr / R2));
 195     D3Ndir.Subtract (
 196       (DNdir.Multiplied ((3.0 * theOffset) * ((D2r/R2) + (Dr*Dr)/R4))));
 197     D3Ndir.Add (Ndir.Multiplied (
 198       (theOffset * (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r))));
 199     D3Ndir.Multiply (theOffset/R);
 200     // V2 = P" (U) :
 201     Standard_Real R4 = R2 * R2;
 202     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
 203     D2Ndir.Subtract (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
 204     D2Ndir.Multiply (theOffset / R);
 205     // V1 = P' (U) :
 206     DNdir.Multiply(R);
 207     DNdir.Subtract (Ndir.Multiplied (Dr/R));
 208     DNdir.Multiply (theOffset/R2);
 209   }
 210   else
 211   {
 212     // Same computation as IICURV in EUCLID-IS because the stability is better.
 213     // V3 = P"' (U) :
 214     D3Ndir.Multiply (theOffset/R);
 215     D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * theOffset * Dr / R3));
 216     D3Ndir.Subtract (DNdir.Multiplied (
 217       ((3.0 * theOffset) * ((D2r/R3) + (Dr*Dr)/R5))) );
 218     D3Ndir.Add (Ndir.Multiplied (
 219       (theOffset * (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r))));
 220     // V2 = P" (U) :
 221     D2Ndir.Multiply (theOffset/R);
 222     D2Ndir.Subtract (DNdir.Multiplied (2.0 * theOffset * Dr / R3));
 223     D2Ndir.Subtract (Ndir.Multiplied (
 224       theOffset * (((3.0 * Dr * Dr) / R5) - (D2r / R3))));
 225     // V1 = P' (U) :
 226     DNdir.Multiply (theOffset/R);
 227     DNdir.Subtract (Ndir.Multiplied (theOffset*Dr/R3));
 228   }
 229
 230   // P(u)
 231   D0(theBasePoint, theBaseD1, theOffset, theIsDirectionChange, theResPoint);
 232   // P'(u)
 233   theResD1 = theBaseD1.Added(gp_Vec2d(DNdir));
 234   // P"(u)
 235   theResD2 = theBaseD2.Added(gp_Vec2d(D2Ndir));
 236   // P"'(u)
 237   if (theIsDirectionChange)
 238     theResD3 = -theBaseD3;
 239   else
 240     theResD3 = theBaseD3;
 241   theResD3.Add(gp_Vec2d(D2Ndir));
 242 }
 243
 244
 245 // ==========  Offset values for 3D curves  ==========
 246
 247 void CSLib_Offset::D0(const gp_Pnt&    theBasePoint,
 248                       const gp_Vec&    theBaseDeriv,
 249                       const gp_Dir&    theOffsetDirection,
 250                       Standard_Real    theOffsetValue,
 251                       Standard_Boolean , // unused
 252                       gp_Pnt&          theResPoint)
 253 {
 254   gp_XYZ Ndir = (theBaseDeriv.XYZ()).Crossed(theOffsetDirection.XYZ());
 255   Standard_Real R = Ndir.Modulus();
 256   if (R <= gp::Resolution())
 257     Standard_NullValue::Raise("CSLib_Offset: Undefined normal vector "
 258           "because tangent vector has zero-magnitude!");
 259
 260   Ndir.Multiply(theOffsetValue / R);
 261   Ndir.Add(theBasePoint.XYZ());
 262   theResPoint.SetXYZ(Ndir);
 263 }
 264
 265 void CSLib_Offset::D1(const gp_Pnt&   theBasePoint,
 266                       const gp_Vec&   theBaseD1,
 267                       const gp_Vec&   theBaseD2,
 268                       const gp_Dir&   theOffsetDirection,
 269                       Standard_Real   theOffsetValue,
 270                       Standard_Boolean , // unused
 271                       gp_Pnt&         theResPoint,
 272                       gp_Vec&         theResDeriv)
 273 {
 274   // P(u) = p(u) + Offset * Ndir / R
 275   // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
 276
 277   // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 278
 279   gp_XYZ Ndir  = (theBaseD1.XYZ()).Crossed(theOffsetDirection.XYZ());
 280   gp_XYZ DNdir = (theBaseD2.XYZ()).Crossed(theOffsetDirection.XYZ());
 281   Standard_Real R2 = Ndir.SquareModulus();
 282   Standard_Real R  = Sqrt (R2);
 283   Standard_Real R3 = R * R2;
 284   Standard_Real Dr = Ndir.Dot (DNdir);
 285   if (R3 <= gp::Resolution()) {
 286     if (R2 <= gp::Resolution())
 287       Standard_NullValue::Raise("CSLib_Offset: Null derivative");
 288     //We try another computation but the stability is not very good.
 289     DNdir.Multiply(R);
 290     DNdir.Subtract(Ndir.Multiplied(Dr / R));
 291     DNdir.Multiply(theOffsetValue / R2);
 292   }
 293   else {
 294     // Same computation as IICURV in EUCLID-IS because the stability is
 295     // better
 296     DNdir.Multiply(theOffsetValue / R);
 297     DNdir.Subtract(Ndir.Multiplied(theOffsetValue * Dr / R3));
 298   }
 299
 300   // P(u)
 301   D0(theBasePoint, theBaseD1, theOffsetDirection, theOffsetValue, Standard_False, theResPoint);
 302   // P'(u)
 303   theResDeriv = theBaseD1.Added(gp_Vec(DNdir));
 304 }
 305
 306 void CSLib_Offset::D2(const gp_Pnt&    theBasePoint,
 307                       const gp_Vec&    theBaseD1,
 308                       const gp_Vec&    theBaseD2,
 309                       const gp_Vec&    theBaseD3,
 310                       const gp_Dir&    theOffsetDirection,
 311                       Standard_Real    theOffsetValue,
 312                       Standard_Boolean theIsDirectionChange,
 313                       gp_Pnt&          theResPoint,
 314                       gp_Vec&          theResD1,
 315                       gp_Vec&          theResD2)
 316 {
 317   // P(u) = p(u) + Offset * Ndir / R
 318   // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
 319
 320   // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 321
 322   // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
 323   //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 324
 325   gp_XYZ Ndir   = (theBaseD1.XYZ()).Crossed(theOffsetDirection.XYZ());
 326   gp_XYZ DNdir  = (theBaseD2.XYZ()).Crossed(theOffsetDirection.XYZ());
 327   gp_XYZ D2Ndir = (theBaseD3.XYZ()).Crossed(theOffsetDirection.XYZ());
 328   Standard_Real R2  = Ndir.SquareModulus();
 329   Standard_Real R   = Sqrt (R2);
 330   Standard_Real R3  = R2 * R;
 331   Standard_Real R4  = R2 * R2;
 332   Standard_Real R5  = R3 * R2;
 333   Standard_Real Dr  = Ndir.Dot (DNdir);
 334   Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
 335
 336   if (R5 <= gp::Resolution()) {
 337     if (R4 <= gp::Resolution())
 338       Standard_NullValue::Raise("CSLib_Offset: Null derivative");
 339     //We try another computation but the stability is not very good
 340     //dixit ISG.
 341     // V2 = P" (U) :
 342     Standard_Real R4 = R2 * R2;
 343     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
 344     D2Ndir.Add (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
 345     D2Ndir.Multiply (theOffsetValue / R);
 346
 347     // V1 = P' (U) :
 348     DNdir.Multiply(R);
 349     DNdir.Subtract (Ndir.Multiplied (Dr/R));
 350     DNdir.Multiply (theOffsetValue/R2);
 351   }
 352   else {
 353     // Same computation as IICURV in EUCLID-IS because the stability is
 354     // better.
 355     // V2 = P" (U) :
 356     D2Ndir.Multiply (theOffsetValue/R);
 357     D2Ndir.Subtract (DNdir.Multiplied (2.0 * theOffsetValue * Dr / R3));
 358     D2Ndir.Add (Ndir.Multiplied (theOffsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))));
 359
 360     // V1 = P' (U) :
 361     DNdir.Multiply (theOffsetValue/R);
 362     DNdir.Subtract (Ndir.Multiplied (theOffsetValue*Dr/R3));
 363   }
 364
 365   // P(u) :
 366   D0(theBasePoint, theBaseD1, theOffsetDirection, theOffsetValue, theIsDirectionChange, theResPoint);
 367   // P'(u) :
 368   theResD1 = theBaseD1.Added(gp_Vec(DNdir));
 369   // P"(u) :
 370   if (theIsDirectionChange)
 371     theResD2 = -theBaseD2;
 372   else
 373     theResD2 = theBaseD2;
 374   theResD2.Add(gp_Vec(D2Ndir));
 375 }
 376
 377 void CSLib_Offset::D3(const gp_Pnt&    theBasePoint,
 378                       const gp_Vec&    theBaseD1,
 379                       const gp_Vec&    theBaseD2,
 380                       const gp_Vec&    theBaseD3,
 381                       const gp_Vec&    theBaseD4,
 382                       const gp_Dir&    theOffsetDirection,
 383                       Standard_Real    theOffsetValue,
 384                       Standard_Boolean theIsDirectionChange,
 385                       gp_Pnt&          theResPoint,
 386                       gp_Vec&          theResD1,
 387                       gp_Vec&          theResD2,
 388                       gp_Vec&          theResD3)
 389 {
 390   // P(u) = p(u) + Offset * Ndir / R
 391   // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
 392
 393   // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 394
 395   // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
 396   //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 397
 398   //P"'(u) = p"'(u) + (Offset / R) * (D3Ndir - (3.0 * Dr/R**2) * D2Ndir -
 399   //         (3.0 * D2r / R2) * DNdir + (3.0 * Dr * Dr / R4) * DNdir -
 400   //         (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir +
 401   //         (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir
 402
 403   gp_XYZ Ndir   = (theBaseD1.XYZ()).Crossed(theOffsetDirection.XYZ());
 404   gp_XYZ DNdir  = (theBaseD2.XYZ()).Crossed(theOffsetDirection.XYZ());
 405   gp_XYZ D2Ndir = (theBaseD3.XYZ()).Crossed(theOffsetDirection.XYZ());
 406   gp_XYZ D3Ndir = (theBaseD4.XYZ()).Crossed(theOffsetDirection.XYZ());
 407   Standard_Real R2  = Ndir.SquareModulus();
 408   Standard_Real R   = Sqrt (R2);
 409   Standard_Real R3  = R2 * R;
 410   Standard_Real R4  = R2 * R2;
 411   Standard_Real R5  = R3 * R2;
 412   Standard_Real R6  = R3 * R3;
 413   Standard_Real R7  = R5 * R2;
 414   Standard_Real Dr  = Ndir.Dot (DNdir);
 415   Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
 416   Standard_Real D3r = Ndir.Dot (D3Ndir) + 3.0 * DNdir.Dot (D2Ndir);
 417   if (R7 <= gp::Resolution()) {
 418     if (R6 <= gp::Resolution())
 419       Standard_NullValue::Raise("CSLib_Offset: Null derivative");
 420     // V3 = P"' (U) :
 421     D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * Dr / R2));
 422     D3Ndir.Subtract (DNdir.Multiplied (3.0 * ((D2r/R2) + (Dr*Dr/R4))));
 423     D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r));
 424     D3Ndir.Multiply (theOffsetValue/R);
 425     // V2 = P" (U) :
 426     Standard_Real R4 = R2 * R2;
 427     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
 428     D2Ndir.Subtract (Ndir.Multiplied ((3.0 * Dr * Dr / R4) - (D2r / R2)));
 429     D2Ndir.Multiply (theOffsetValue / R);
 430     // V1 = P' (U) :
 431     DNdir.Multiply(R);
 432     DNdir.Subtract (Ndir.Multiplied (Dr/R));
 433     DNdir.Multiply (theOffsetValue/R2);
 434   }
 435   else {
 436     // V3 = P"' (U) :
 437     D3Ndir.Divide (R);
 438     D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * Dr / R3));
 439     D3Ndir.Subtract (DNdir.Multiplied ((3.0 * ((D2r/R3) + (Dr*Dr)/R5))));
 440     D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r));
 441     D3Ndir.Multiply (theOffsetValue);
 442     // V2 = P" (U) :
 443     D2Ndir.Divide (R);
 444     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R3));
 445     D2Ndir.Subtract (Ndir.Multiplied ((3.0 * Dr * Dr / R5) - (D2r / R3)));
 446     D2Ndir.Multiply (theOffsetValue);
 447     // V1 = P' (U) :
 448     DNdir.Multiply (theOffsetValue/R);
 449     DNdir.Subtract (Ndir.Multiplied (theOffsetValue*Dr/R3));
 450   }
 451
 452   // P(u)
 453   D0(theBasePoint, theBaseD1, theOffsetDirection, theOffsetValue, theIsDirectionChange, theResPoint);
 454   // P'(u)
 455   theResD1 = theBaseD1.Added(gp_Vec(DNdir));
 456   // P"(u)
 457   theResD2 = theBaseD2.Added(gp_Vec(D2Ndir));
 458   // P"'(u)
 459   if (theIsDirectionChange)
 460     theResD3 = -theBaseD3;
 461   else
 462     theResD3 = theBaseD3;
 463   theResD3.Add(gp_Vec(D2Ndir));
 464 }
 465