0029948: Uninitialized variable in GeomEvaluator_OffsetSurface::CalculateD0(...)...
[occt.git] / src / GeomEvaluator / GeomEvaluator_OffsetSurface.cxx
index acc7e30..6943c92 100644 (file)
 #include <gp_Vec.hxx>
 #include <Standard_RangeError.hxx>
 
+IMPLEMENT_STANDARD_RTTIEXT(GeomEvaluator_OffsetSurface,GeomEvaluator_Surface)
 
+namespace {
+
+// tolerance for considering derivative to be null
+const Standard_Real the_D1MagTol = 1.e-9;
+
+// If calculation of normal fails, try shifting the point towards the center 
+// of the parametric space of the surface, in the hope that derivatives 
+// are better defined there.
+//
+// This shift is iterative, starting with Precision::PConfusion()
+// and increasing by multiple of 2 on each step.
+//
+// NB: temporarily this is made as static function and not class method, 
+// hence code duplications
+static Standard_Boolean shiftPoint (const Standard_Real theUStart, const Standard_Real theVStart,
+                                    Standard_Real& theU, Standard_Real& theV, 
+                                    const Handle(Geom_Surface)& theSurf,
+                                    const Handle(GeomAdaptor_HSurface)& theAdaptor,
+                                    const gp_Vec& theD1U, const gp_Vec& theD1V)
+{
+  // Get parametric bounds and closure status
+  Standard_Real aUMin, aUMax, aVMin, aVMax;
+  Standard_Boolean isUPeriodic, isVPeriodic;
+  if (! theSurf.IsNull())
+  {
+    theSurf->Bounds (aUMin, aUMax, aVMin, aVMax);
+    isUPeriodic = theSurf->IsUPeriodic();
+    isVPeriodic = theSurf->IsVPeriodic();
+  }
+  else
+  {
+    aUMin = theAdaptor->FirstUParameter();
+    aUMax = theAdaptor->LastUParameter();
+    aVMin = theAdaptor->FirstVParameter();
+    aVMax = theAdaptor->LastVParameter();
+    isUPeriodic = theAdaptor->IsUPeriodic();
+    isVPeriodic = theAdaptor->IsVPeriodic();
+  }
+
+  // check if either U or V is singular (normally one of them is)
+  Standard_Boolean isUSingular = (theD1U.SquareMagnitude() < the_D1MagTol * the_D1MagTol);
+  Standard_Boolean isVSingular = (theD1V.SquareMagnitude() < the_D1MagTol * the_D1MagTol);
+
+  // compute vector to shift from start point to center of the surface;
+  // if surface is periodic or singular in some direction, take shift in that direction zero
+  Standard_Real aDirU = (isUPeriodic || (isUSingular && ! isVSingular) ? 0. : 0.5 * (aUMin + aUMax) - theUStart);
+  Standard_Real aDirV = (isVPeriodic || (isVSingular && ! isUSingular) ? 0. : 0.5 * (aVMin + aVMax) - theVStart);
+  Standard_Real aDist = Sqrt (aDirU * aDirU + aDirV * aDirV);
+
+  // shift current point from its current position towards center, by value of twice
+  // current distance from it to start (but not less than Precision::PConfusion());
+  // fail if center is overpassed.
+  Standard_Real aDU = theU - theUStart;
+  Standard_Real aDV = theV - theVStart;
+  Standard_Real aStep = Max (2. * Sqrt (aDU * aDU + aDV * aDV), Precision::PConfusion());
+  if (aStep >= aDist)
+  {
+    return Standard_False;
+  }
+
+  aStep /= aDist;
+  theU += aDirU * aStep;
+  theV += aDirV * aStep;
+
+//  std::cout << "Normal calculation failed at (" << theUStart << ", " << theVStart << "), shifting to (" << theU << ", " << theV << ")" << std::endl;
+
+  return Standard_True;
+}
+
+// auxiliary function
 template<class SurfOrAdapt>
 static void derivatives(Standard_Integer theMaxOrder,
                         Standard_Integer theMinOrder,
@@ -37,9 +108,103 @@ static void derivatives(Standard_Integer theMaxOrder,
                         const Standard_Boolean theAlongV,
                         const Handle(Geom_BSplineSurface)& theL,
                         TColgp_Array2OfVec& theDerNUV,
-                        TColgp_Array2OfVec& theDerSurf);
+                        TColgp_Array2OfVec& theDerSurf)
+{
+  Standard_Integer i, j;
+  gp_Pnt P;
+  gp_Vec DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3UUV, DL3UVV, DL3V;
 
+  if (theAlongU || theAlongV)
+  {
+    theMaxOrder = 0;
+    TColgp_Array2OfVec DerSurfL(0, theMaxOrder + theNU + 1, 0, theMaxOrder + theNV + 1);
+    switch (theMinOrder)
+    {
+    case 1:
+      theL->D1(theU, theV, P, DL1U, DL1V);
+      DerSurfL.SetValue(1, 0, DL1U);
+      DerSurfL.SetValue(0, 1, DL1V);
+      break;
+    case 2:
+      theL->D2(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV);
+      DerSurfL.SetValue(1, 0, DL1U);
+      DerSurfL.SetValue(0, 1, DL1V);
+      DerSurfL.SetValue(1, 1, DL2UV);
+      DerSurfL.SetValue(2, 0, DL2U);
+      DerSurfL.SetValue(0, 2, DL2V);
+      break;
+    case 3:
+      theL->D3(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3V, DL3UUV, DL3UVV);
+      DerSurfL.SetValue(1, 0, DL1U);
+      DerSurfL.SetValue(0, 1, DL1V);
+      DerSurfL.SetValue(1, 1, DL2UV);
+      DerSurfL.SetValue(2, 0, DL2U);
+      DerSurfL.SetValue(0, 2, DL2V);
+      DerSurfL.SetValue(3, 0, DL3U);
+      DerSurfL.SetValue(2, 1, DL3UUV);
+      DerSurfL.SetValue(1, 2, DL3UVV);
+      DerSurfL.SetValue(0, 3, DL3V);
+      break;
+    default:
+      break;
+    }
+
+    if (theNU <= theNV)
+    {
+      for (i = 0; i <= theMaxOrder + 1 + theNU; i++)
+        for (j = i; j <= theMaxOrder + theNV + 1; j++)
+          if (i + j > theMinOrder)
+          {
+            DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
+            theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
+            if (i != j && j <= theNU + 1)
+            {
+              theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
+              DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
+            }
+          }
+    }
+    else
+    {
+      for (j = 0; j <= theMaxOrder + 1 + theNV; j++)
+        for (i = j; i <= theMaxOrder + theNU + 1; i++)
+          if (i + j > theMinOrder)
+          {
+            DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
+            theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
+            if (i != j && i <= theNV + 1)
+            {
+              theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
+              DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
+            }
+          }
+    }
+    for (i = 0; i <= theMaxOrder + theNU; i++)
+      for (j = 0; j <= theMaxOrder + theNV; j++)
+      {
+        if (theAlongU)
+          theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, DerSurfL, theDerSurf));
+        if (theAlongV)
+          theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf, DerSurfL));
+      }
+  }
+  else
+  {
+    for (i = 0; i <= theMaxOrder + theNU+ 1; i++)
+      for (j = i; j <= theMaxOrder + theNV + 1; j++)
+        if (i + j > theMinOrder)
+        {
+          theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
+          if (i != j)
+            theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
+        }
+    for (i = 0; i <= theMaxOrder + theNU; i++)
+      for (j = 0; j <= theMaxOrder + theNV; j++)
+        theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf));
+  }
+}
 
+} // end of anonymous namespace
 
 GeomEvaluator_OffsetSurface::GeomEvaluator_OffsetSurface(
         const Handle(Geom_Surface)& theBase,
@@ -65,39 +230,58 @@ GeomEvaluator_OffsetSurface::GeomEvaluator_OffsetSurface(
         const Handle(Geom_OsculatingSurface)& theOscSurf)
   : GeomEvaluator_Surface(),
     myBaseAdaptor(theBase),
-    myOffset(theOffset)
+    myOffset(theOffset),
+    myOscSurf(theOscSurf)
 {
-  if (theOscSurf.IsNull())
-    return; // osculating surface already exists
-
-  // Create osculating surface for B-spline and Besier surfaces only
-  Handle(Geom_Surface) aBSurf;
-  if (myBaseAdaptor->GetType() == GeomAbs_BSplineSurface)
-    aBSurf = myBaseAdaptor->BSpline();
-  else if (myBaseAdaptor->GetType() == GeomAbs_BezierSurface)
-    aBSurf = myBaseAdaptor->Bezier();
-  if (!aBSurf.IsNull())
-  {
-    myBaseSurf = aBSurf;
-    myOscSurf = new Geom_OsculatingSurface(aBSurf, Precision::Confusion());
-  }
 }
 
 void GeomEvaluator_OffsetSurface::D0(
     const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue) const
 {
-  gp_Vec aD1U, aD1V;
-  BaseD1(theU, theV, theValue, aD1U, aD1V);
-  CalculateD0(theU, theV, theValue, aD1U, aD1V);
+  Standard_Real aU = theU, aV = theV;
+  for (;;)
+  {
+    gp_Vec aD1U, aD1V;
+    BaseD1 (aU, aV, theValue, aD1U, aD1V);
+    try 
+    {
+      CalculateD0(aU, aV, theValue, aD1U, aD1V);
+      break;
+    } 
+    catch (Geom_UndefinedValue&)
+    {
+      // if failed at parametric boundary, try taking derivative at shifted point
+      if (! shiftPoint (theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, aD1U, aD1V))
+      {
+        throw;
+      }
+    }
+  }
 }
 
 void GeomEvaluator_OffsetSurface::D1(
     const Standard_Real theU, const Standard_Real theV,
     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V) const
 {
-  gp_Vec aD2U, aD2V, aD2UV;
-  BaseD2(theU, theV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
-  CalculateD1(theU, theV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
+  Standard_Real aU = theU, aV = theV;
+  for (;;)
+  {
+    gp_Vec aD2U, aD2V, aD2UV;
+    BaseD2 (aU, aV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
+    try 
+    {
+      CalculateD1(aU, aV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
+      break;
+    } 
+    catch (Geom_UndefinedValue&)
+    {
+      // if failed at parametric boundary, try taking derivative at shifted point
+      if (! shiftPoint (theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, theD1U, theD1V))
+      {
+        throw;
+      }
+    }
+  }
 }
 
 void GeomEvaluator_OffsetSurface::D2(
@@ -105,11 +289,27 @@ void GeomEvaluator_OffsetSurface::D2(
     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV) const
 {
-  gp_Vec aD3U, aD3V, aD3UUV, aD3UVV;
-  BaseD3(theU, theV, theValue, theD1U, theD1V,
-         theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
-  CalculateD2(theU, theV, theValue, theD1U, theD1V,
-              theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
+  Standard_Real aU = theU, aV = theV;
+  for (;;)
+  {
+    gp_Vec aD3U, aD3V, aD3UUV, aD3UVV;
+    BaseD3 (aU, aV, theValue, theD1U, theD1V,
+            theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
+    try
+    {
+      CalculateD2 (aU, aV, theValue, theD1U, theD1V,
+                   theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
+      break;
+    } 
+    catch (Geom_UndefinedValue&)
+    {
+      // if failed at parametric boundary, try taking derivative at shifted point
+      if (! shiftPoint (theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, theD1U, theD1V))
+      {
+        throw;
+      }
+    }
+  }
 }
 
 void GeomEvaluator_OffsetSurface::D3(
@@ -118,10 +318,26 @@ void GeomEvaluator_OffsetSurface::D3(
     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
     gp_Vec& theD3U, gp_Vec& theD3V, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const
 {
-  BaseD3(theU, theV, theValue, theD1U, theD1V,
-         theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
-  CalculateD3(theU, theV, theValue, theD1U, theD1V,
-              theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
+  Standard_Real aU = theU, aV = theV;
+  for (;;)
+  {
+    BaseD3 (aU, aV, theValue, theD1U, theD1V,
+            theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
+    try
+    {
+      CalculateD3 (aU, aV, theValue, theD1U, theD1V,
+                   theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
+      break;
+    } 
+    catch (Geom_UndefinedValue&)
+    {
+      // if failed at parametric boundary, try taking derivative at shifted point
+      if (! shiftPoint (theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, theD1U, theD1V))
+      {
+        throw;
+      }
+    }
+  }
 }
 
 gp_Vec GeomEvaluator_OffsetSurface::DN(
@@ -133,10 +349,25 @@ gp_Vec GeomEvaluator_OffsetSurface::DN(
   Standard_RangeError_Raise_if(theDerU + theDerV < 1,
       "GeomEvaluator_OffsetSurface::DN(): theDerU + theDerV < 1");
 
-  gp_Pnt aP;
-  gp_Vec aD1U, aD1V;
-  BaseD1(theU, theV, aP, aD1U, aD1V);
-  return CalculateDN(theU, theV, theDerU, theDerV, aD1U, aD1V);
+  Standard_Real aU = theU, aV = theV;
+  for (;;)
+  {
+    gp_Pnt aP;
+    gp_Vec aD1U, aD1V;
+    BaseD1 (aU, aV, aP, aD1U, aD1V);
+    try
+    {
+      return CalculateDN (aU, aV, theDerU, theDerV, aD1U, aD1V);
+    } 
+    catch (Geom_UndefinedValue&)
+    {
+      // if failed at parametric boundary, try taking derivative at shifted point
+      if (! shiftPoint (theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, aD1U, aD1V))
+      {
+        throw;
+      }
+    }
+  }
 }
 
 
@@ -188,8 +419,6 @@ void GeomEvaluator_OffsetSurface::CalculateD0(
     gp_Pnt& theValue,
     const gp_Vec& theD1U, const gp_Vec& theD1V) const
 {
-  const Standard_Real MagTol = 1.e-9;
-
   // Normalize derivatives before normal calculation because it gives more stable result.
   // There will be normalized only derivatives greater than 1.0 to avoid differences in last significant digit
   gp_Vec aD1U(theD1U);
@@ -202,7 +431,7 @@ void GeomEvaluator_OffsetSurface::CalculateD0(
     aD1V /= Sqrt(aD1VNorm2);
 
   gp_Vec aNorm = aD1U.Crossed(aD1V);
-  if (aNorm.SquareMagnitude() > MagTol * MagTol)
+  if (aNorm.SquareMagnitude() > the_D1MagTol * the_D1MagTol)
   {
     // Non singular case. Simple computations.
     aNorm.Normalize();
@@ -237,20 +466,20 @@ void GeomEvaluator_OffsetSurface::CalculateD0(
       derivatives(MaxOrder, 1, theU, theV, myBaseAdaptor, 0, 0, AlongU, AlongV, L, DerNUV, DerSurf);
 
     gp_Dir Normal;
-    CSLib_NormalStatus NStatus;
-    CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
+    CSLib_NormalStatus NStatus = CSLib_Singular;
+    CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
                   NStatus, Normal, OrderU, OrderV);
     if (NStatus == CSLib_InfinityOfSolutions)
     {
       // Replace zero derivative and try to calculate normal
       gp_Vec aNewDU = theD1U;
       gp_Vec aNewDV = theD1V;
-      if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, MagTol * MagTol))
-        CSLib::Normal(aNewDU, aNewDV, MagTol, NStatus, Normal);
+      if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, the_D1MagTol * the_D1MagTol))
+        CSLib::Normal(aNewDU, aNewDV, the_D1MagTol, NStatus, Normal);
     }
 
     if (NStatus != CSLib_Defined)
-      Geom_UndefinedValue::Raise(
+      throw Geom_UndefinedValue(
           "GeomEvaluator_OffsetSurface::CalculateD0(): Unable to calculate normal");
 
     theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
@@ -262,19 +491,11 @@ void GeomEvaluator_OffsetSurface::CalculateD1(
     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
     const gp_Vec& theD2U, const gp_Vec& theD2V, const gp_Vec& theD2UV) const
 {
-  const Standard_Real MagTol = 1.e-9;
-
   // Check offset side.
   Handle(Geom_BSplineSurface) L;
   Standard_Boolean isOpposite = Standard_False;
   Standard_Boolean AlongU = Standard_False;
   Standard_Boolean AlongV = Standard_False;
-  if (!myOscSurf.IsNull())
-  {
-    AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
-    AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
-  }
-  const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
 
   // Normalize derivatives before normal calculation because it gives more stable result.
   // There will be normalized only derivatives greater than 1.0 to avoid differences in last significant digit
@@ -287,48 +508,58 @@ void GeomEvaluator_OffsetSurface::CalculateD1(
   if (aD1VNorm2 > 1.0)
     aD1V /= Sqrt(aD1VNorm2);
 
+  Standard_Boolean isSingular = Standard_False;
   Standard_Integer MaxOrder = 0;
   gp_Vec aNorm = aD1U.Crossed(aD1V);
-  if (aNorm.SquareMagnitude() > MagTol * MagTol)
+  if (aNorm.SquareMagnitude() <= the_D1MagTol * the_D1MagTol)
   {
-    if (!AlongV && !AlongU)
+    if (!myOscSurf.IsNull())
     {
-      // AlongU or AlongV leads to more complex D1 computation
-      // Try to compute D0 and D1 much simpler
-      aNorm.Normalize();
-      theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * aNorm.XYZ());
-
-      gp_Vec aN0(aNorm.XYZ()), aN1U, aN1V;
-      Standard_Real aScale = (theD1U^theD1V).Dot(aN0);
-      aN1U.SetX(theD2U.Y() * theD1V.Z() + theD1U.Y() * theD2UV.Z()
-        - theD2U.Z() * theD1V.Y() - theD1U.Z() * theD2UV.Y());
-      aN1U.SetY((theD2U.X() * theD1V.Z() + theD1U.X() * theD2UV.Z()
-        - theD2U.Z() * theD1V.X() - theD1U.Z() * theD2UV.X()) * -1.0);
-      aN1U.SetZ(theD2U.X() * theD1V.Y() + theD1U.X() * theD2UV.Y()
-        - theD2U.Y() * theD1V.X() - theD1U.Y() * theD2UV.X());
-      Standard_Real aScaleU = aN1U.Dot(aN0);
-      aN1U.Subtract(aScaleU * aN0);
-      aN1U /= aScale;
-
-      aN1V.SetX(theD2UV.Y() * theD1V.Z() + theD2V.Z() * theD1U.Y()
-        - theD2UV.Z() * theD1V.Y() - theD2V.Y() * theD1U.Z());
-      aN1V.SetY((theD2UV.X() * theD1V.Z() + theD2V.Z() * theD1U.X()
-        - theD2UV.Z() * theD1V.X() - theD2V.X() * theD1U.Z()) * -1.0);
-      aN1V.SetZ(theD2UV.X() * theD1V.Y() + theD2V.Y() * theD1U.X()
-        - theD2UV.Y() * theD1V.X() - theD2V.X() * theD1U.Y());
-      Standard_Real aScaleV = aN1V.Dot(aN0);
-      aN1V.Subtract(aScaleV * aN0);
-      aN1V /= aScale;
-
-      theD1U += myOffset * aSign * aN1U;
-      theD1V += myOffset * aSign * aN1V;
-
-      return;
+      AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
+      AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
     }
-  }
-  else
+
     MaxOrder = 3;
+    isSingular = Standard_True;
+  }
 
+  const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
+  
+  if (!isSingular)
+  {
+    // AlongU or AlongV leads to more complex D1 computation
+    // Try to compute D0 and D1 much simpler
+    aNorm.Normalize();
+    theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * aNorm.XYZ());
+
+    gp_Vec aN0(aNorm.XYZ()), aN1U, aN1V;
+    Standard_Real aScale = (theD1U^theD1V).Dot(aN0);
+    aN1U.SetX(theD2U.Y() * theD1V.Z() + theD1U.Y() * theD2UV.Z()
+      - theD2U.Z() * theD1V.Y() - theD1U.Z() * theD2UV.Y());
+    aN1U.SetY((theD2U.X() * theD1V.Z() + theD1U.X() * theD2UV.Z()
+      - theD2U.Z() * theD1V.X() - theD1U.Z() * theD2UV.X()) * -1.0);
+    aN1U.SetZ(theD2U.X() * theD1V.Y() + theD1U.X() * theD2UV.Y()
+      - theD2U.Y() * theD1V.X() - theD1U.Y() * theD2UV.X());
+    Standard_Real aScaleU = aN1U.Dot(aN0);
+    aN1U.Subtract(aScaleU * aN0);
+    aN1U /= aScale;
+
+    aN1V.SetX(theD2UV.Y() * theD1V.Z() + theD2V.Z() * theD1U.Y()
+      - theD2UV.Z() * theD1V.Y() - theD2V.Y() * theD1U.Z());
+    aN1V.SetY((theD2UV.X() * theD1V.Z() + theD2V.Z() * theD1U.X()
+      - theD2UV.Z() * theD1V.X() - theD2V.X() * theD1U.Z()) * -1.0);
+    aN1V.SetZ(theD2UV.X() * theD1V.Y() + theD2V.Y() * theD1U.X()
+      - theD2UV.Y() * theD1V.X() - theD2V.X() * theD1U.Y());
+    Standard_Real aScaleV = aN1V.Dot(aN0);
+    aN1V.Subtract(aScaleV * aN0);
+    aN1V /= aScale;
+
+    theD1U += myOffset * aSign * aN1U;
+    theD1V += myOffset * aSign * aN1V;
+
+    return;
+  }
+  
   Standard_Integer OrderU, OrderV;
   TColgp_Array2OfVec DerNUV(0, MaxOrder + 1, 0, MaxOrder + 1);
   TColgp_Array2OfVec DerSurf(0, MaxOrder + 2, 0, MaxOrder + 2);
@@ -347,13 +578,13 @@ void GeomEvaluator_OffsetSurface::CalculateD1(
 
   gp_Dir Normal;
   CSLib_NormalStatus NStatus;
-  CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax, NStatus, Normal, OrderU, OrderV);
+  CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax, NStatus, Normal, OrderU, OrderV);
   if (NStatus == CSLib_InfinityOfSolutions)
   {
     gp_Vec aNewDU = theD1U;
     gp_Vec aNewDV = theD1V;
     // Replace zero derivative and try to calculate normal
-    if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, MagTol * MagTol))
+    if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, the_D1MagTol * the_D1MagTol))
     {
       DerSurf.SetValue(1, 0, aNewDU);
       DerSurf.SetValue(0, 1, aNewDV);
@@ -361,13 +592,13 @@ void GeomEvaluator_OffsetSurface::CalculateD1(
         derivatives(MaxOrder, 2, theU, theV, myBaseSurf, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
       else
         derivatives(MaxOrder, 2, theU, theV, myBaseAdaptor, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
-      CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
+      CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
                     NStatus, Normal, OrderU, OrderV);
     }
   }
 
   if (NStatus != CSLib_Defined)
-    Geom_UndefinedValue::Raise(
+    throw Geom_UndefinedValue(
         "GeomEvaluator_OffsetSurface::CalculateD1(): Unable to calculate normal");
 
   theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
@@ -382,11 +613,9 @@ void GeomEvaluator_OffsetSurface::CalculateD2(
     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
     const gp_Vec& theD3U, const gp_Vec& theD3V, const gp_Vec& theD3UUV, const gp_Vec& theD3UVV) const
 {
-  const Standard_Real MagTol = 1.e-9;
-
   gp_Dir Normal;
   CSLib_NormalStatus NStatus;
-  CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
+  CSLib::Normal(theD1U, theD1V, the_D1MagTol, NStatus, Normal);
 
   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
   Standard_Integer OrderU, OrderV;
@@ -411,7 +640,7 @@ void GeomEvaluator_OffsetSurface::CalculateD2(
   Standard_Boolean isOpposite = Standard_False;
   Standard_Boolean AlongU = Standard_False;
   Standard_Boolean AlongV = Standard_False;
-  if (!myOscSurf.IsNull())
+  if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull())
   {
     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
@@ -423,10 +652,10 @@ void GeomEvaluator_OffsetSurface::CalculateD2(
   else
     derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 2, 2, AlongU, AlongV, L, DerNUV, DerSurf);
 
-  CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
+  CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
                 NStatus, Normal, OrderU, OrderV);
   if (NStatus != CSLib_Defined)
-    Geom_UndefinedValue::Raise(
+    throw Geom_UndefinedValue(
         "GeomEvaluator_OffsetSurface::CalculateD2(): Unable to calculate normal");
 
   theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
@@ -458,11 +687,9 @@ void GeomEvaluator_OffsetSurface::CalculateD3(
     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
     gp_Vec& theD3U, gp_Vec& theD3V, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const
 {
-  const Standard_Real MagTol = 1.e-9;
-
   gp_Dir Normal;
   CSLib_NormalStatus NStatus;
-  CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
+  CSLib::Normal(theD1U, theD1V, the_D1MagTol, NStatus, Normal);
   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
   Standard_Integer OrderU, OrderV;
   TColgp_Array2OfVec DerNUV(0, MaxOrder + 3, 0, MaxOrder + 3);
@@ -486,7 +713,7 @@ void GeomEvaluator_OffsetSurface::CalculateD3(
   Standard_Boolean isOpposite = Standard_False;
   Standard_Boolean AlongU = Standard_False;
   Standard_Boolean AlongV = Standard_False;
-  if (!myOscSurf.IsNull())
+  if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull())
   {
     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
@@ -498,10 +725,10 @@ void GeomEvaluator_OffsetSurface::CalculateD3(
   else
     derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 3, 3, AlongU, AlongV, L, DerNUV, DerSurf);
 
-  CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
+  CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
                 NStatus, Normal, OrderU, OrderV);
   if (NStatus != CSLib_Defined)
-    Geom_UndefinedValue::Raise(
+    throw Geom_UndefinedValue(
         "GeomEvaluator_OffsetSurface::CalculateD3(): Unable to calculate normal");
 
   theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
@@ -544,11 +771,9 @@ gp_Vec GeomEvaluator_OffsetSurface::CalculateDN(
     const Standard_Integer theNu, const Standard_Integer theNv,
     const gp_Vec& theD1U, const gp_Vec& theD1V) const
 {
-  const Standard_Real MagTol = 1.e-9;
-
   gp_Dir Normal;
   CSLib_NormalStatus NStatus;
-  CSLib::Normal(theD1U, theD1V, MagTol, NStatus, Normal);
+  CSLib::Normal(theD1U, theD1V, the_D1MagTol, NStatus, Normal);
   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
   Standard_Integer OrderU, OrderV;
   TColgp_Array2OfVec DerNUV(0, MaxOrder + theNu, 0, MaxOrder + theNu);
@@ -566,7 +791,7 @@ gp_Vec GeomEvaluator_OffsetSurface::CalculateDN(
   Standard_Boolean isOpposite = Standard_False;
   Standard_Boolean AlongU = Standard_False;
   Standard_Boolean AlongV = Standard_False;
-  if (!myOscSurf.IsNull())
+  if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull())
   {
     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
@@ -578,10 +803,10 @@ gp_Vec GeomEvaluator_OffsetSurface::CalculateDN(
   else
     derivatives(MaxOrder, 1, theU, theV, myBaseAdaptor, theNu, theNv, AlongU, AlongV, L, DerNUV, DerSurf);
 
-  CSLib::Normal(MaxOrder, DerNUV, MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
+  CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
                 NStatus, Normal, OrderU, OrderV);
   if (NStatus != CSLib_Defined)
-    Geom_UndefinedValue::Raise(
+    throw Geom_UndefinedValue(
         "GeomEvaluator_OffsetSurface::CalculateDN(): Unable to calculate normal");
 
   gp_Vec D;
@@ -680,117 +905,3 @@ Standard_Boolean GeomEvaluator_OffsetSurface::ReplaceDerivative(
   }
   return isReplaced;
 }
-
-
-
-
-
-// =====================   Auxiliary functions   ===================================
-template<class SurfOrAdapt>
-void derivatives(Standard_Integer theMaxOrder,
-                 Standard_Integer theMinOrder,
-                 const Standard_Real theU,
-                 const Standard_Real theV,
-                 const SurfOrAdapt& theBasisSurf,
-                 const Standard_Integer theNU,
-                 const Standard_Integer theNV,
-                 const Standard_Boolean theAlongU,
-                 const Standard_Boolean theAlongV,
-                 const Handle(Geom_BSplineSurface)& theL,
-                 TColgp_Array2OfVec& theDerNUV,
-                 TColgp_Array2OfVec& theDerSurf)
-{
-  Standard_Integer i, j;
-  gp_Pnt P;
-  gp_Vec DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3UUV, DL3UVV, DL3V;
-
-  if (theAlongU || theAlongV)
-  {
-    theMaxOrder = 0;
-    TColgp_Array2OfVec DerSurfL(0, theMaxOrder + theNU + 1, 0, theMaxOrder + theNV + 1);
-    switch (theMinOrder)
-    {
-    case 1:
-      theL->D1(theU, theV, P, DL1U, DL1V);
-      DerSurfL.SetValue(1, 0, DL1U);
-      DerSurfL.SetValue(0, 1, DL1V);
-      break;
-    case 2:
-      theL->D2(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV);
-      DerSurfL.SetValue(1, 0, DL1U);
-      DerSurfL.SetValue(0, 1, DL1V);
-      DerSurfL.SetValue(1, 1, DL2UV);
-      DerSurfL.SetValue(2, 0, DL2U);
-      DerSurfL.SetValue(0, 2, DL2V);
-      break;
-    case 3:
-      theL->D3(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3V, DL3UUV, DL3UVV);
-      DerSurfL.SetValue(1, 0, DL1U);
-      DerSurfL.SetValue(0, 1, DL1V);
-      DerSurfL.SetValue(1, 1, DL2UV);
-      DerSurfL.SetValue(2, 0, DL2U);
-      DerSurfL.SetValue(0, 2, DL2V);
-      DerSurfL.SetValue(3, 0, DL3U);
-      DerSurfL.SetValue(2, 1, DL3UUV);
-      DerSurfL.SetValue(1, 2, DL3UVV);
-      DerSurfL.SetValue(0, 3, DL3V);
-      break;
-    default:
-      break;
-    }
-
-    if (theNU <= theNV)
-    {
-      for (i = 0; i <= theMaxOrder + 1 + theNU; i++)
-        for (j = i; j <= theMaxOrder + theNV + 1; j++)
-          if (i + j > theMinOrder)
-          {
-            DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
-            theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
-            if (i != j && j <= theNU + 1)
-            {
-              theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
-              DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
-            }
-          }
-    }
-    else
-    {
-      for (j = 0; j <= theMaxOrder + 1 + theNV; j++)
-        for (i = j; i <= theMaxOrder + theNU + 1; i++)
-          if (i + j > theMinOrder)
-          {
-            DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
-            theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
-            if (i != j && i <= theNV + 1)
-            {
-              theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
-              DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
-            }
-          }
-    }
-    for (i = 0; i <= theMaxOrder + theNU; i++)
-      for (j = 0; j <= theMaxOrder + theNV; j++)
-      {
-        if (theAlongU)
-          theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, DerSurfL, theDerSurf));
-        if (theAlongV)
-          theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf, DerSurfL));
-      }
-  }
-  else
-  {
-    for (i = 0; i <= theMaxOrder + theNU+ 1; i++)
-      for (j = i; j <= theMaxOrder + theNV + 1; j++)
-        if (i + j > theMinOrder)
-        {
-          theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
-          if (i != j)
-            theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
-        }
-    for (i = 0; i <= theMaxOrder + theNU; i++)
-      for (j = 0; j <= theMaxOrder + theNV; j++)
-        theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf));
-  }
-}
-