0031303: Different calculation of offset direction in Adaptor2d_OffsetCurve and Geom2...
[occt.git] / src / Geom2dEvaluator / Geom2dEvaluator_OffsetCurve.cxx
index a899808..56864c1 100644 (file)
@@ -13,7 +13,7 @@
 // commercial license or contractual agreement.
 
 #include <Geom2dEvaluator_OffsetCurve.hxx>
-
+#include <Geom2dEvaluator.hxx>
 #include <Geom2dAdaptor_HCurve.hxx>
 #include <Standard_NullValue.hxx>
 
@@ -43,7 +43,7 @@ void Geom2dEvaluator_OffsetCurve::D0(const Standard_Real theU,
 {
   gp_Vec2d aD1;
   BaseD1(theU, theValue, aD1);
-  CalculateD0(theValue, aD1);
+  Geom2dEvaluator::CalculateD0(theValue, aD1, myOffset);
 }
 
 void Geom2dEvaluator_OffsetCurve::D1(const Standard_Real theU,
@@ -52,7 +52,7 @@ void Geom2dEvaluator_OffsetCurve::D1(const Standard_Real theU,
 {
   gp_Vec2d aD2;
   BaseD2(theU, theValue, theD1, aD2);
-  CalculateD1(theValue, theD1, aD2);
+  Geom2dEvaluator::CalculateD1(theValue, theD1, aD2, myOffset);
 }
 
 void Geom2dEvaluator_OffsetCurve::D2(const Standard_Real theU,
@@ -70,7 +70,7 @@ void Geom2dEvaluator_OffsetCurve::D2(const Standard_Real theU,
     isDirectionChange = AdjustDerivative(3, theU, theD1, theD2, aD3, aDummyD4);
   }
 
-  CalculateD2(theValue, theD1, theD2, aD3, isDirectionChange);
+  Geom2dEvaluator::CalculateD2(theValue, theD1, theD2, aD3, isDirectionChange, myOffset);
 }
 
 void Geom2dEvaluator_OffsetCurve::D3(const Standard_Real theU,
@@ -86,7 +86,7 @@ void Geom2dEvaluator_OffsetCurve::D3(const Standard_Real theU,
   if (theD1.SquareMagnitude() <= gp::Resolution())
     isDirectionChange = AdjustDerivative(4, theU, theD1, theD2, theD3, aD4);
 
-  CalculateD3(theValue, theD1, theD2, theD3, aD4, isDirectionChange);
+  Geom2dEvaluator::CalculateD3(theValue, theD1, theD2, theD3, aD4, isDirectionChange, myOffset);
 }
 
 gp_Vec2d Geom2dEvaluator_OffsetCurve::DN(const Standard_Real theU,
@@ -184,207 +184,6 @@ gp_Vec2d Geom2dEvaluator_OffsetCurve::BaseDN(const Standard_Real theU,
 }
 
 
-void Geom2dEvaluator_OffsetCurve::CalculateD0(      gp_Pnt2d& theValue,
-                                              const gp_Vec2d& theD1) const
-{
-  if (theD1.SquareMagnitude() <= gp::Resolution())
-    throw Standard_NullValue("Geom2dEvaluator_OffsetCurve: Undefined normal vector "
-                             "because tangent vector has zero-magnitude!");
-
-  gp_Dir2d aNormal(theD1.Y(), -theD1.X());
-  theValue.ChangeCoord().Add(aNormal.XY() * myOffset);
-}
-
-void Geom2dEvaluator_OffsetCurve::CalculateD1(      gp_Pnt2d& theValue,
-                                                    gp_Vec2d& theD1,
-                                              const gp_Vec2d& theD2) const
-{
-  // P(u) = p(u) + Offset * Ndir / R
-  // with R = || p' ^ Z|| and Ndir = P' ^ Z
-
-  // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
-
-  gp_XY Ndir(theD1.Y(), -theD1.X());
-  gp_XY DNdir(theD2.Y(), -theD2.X());
-  Standard_Real R2 = Ndir.SquareModulus();
-  Standard_Real R = Sqrt(R2);
-  Standard_Real R3 = R * R2;
-  Standard_Real Dr = Ndir.Dot(DNdir);
-  if (R3 <= gp::Resolution())
-  {
-    if (R2 <= gp::Resolution())
-      throw Standard_NullValue("Geom2dEvaluator_OffsetCurve: Null derivative");
-    //We try another computation but the stability is not very good.
-    DNdir.Multiply(R);
-    DNdir.Subtract(Ndir.Multiplied(Dr / R));
-    DNdir.Multiply(myOffset / R2);
-  }
-  else
-  {
-    // Same computation as IICURV in EUCLID-IS because the stability is better
-    DNdir.Multiply(myOffset / R);
-    DNdir.Subtract(Ndir.Multiplied(myOffset * Dr / R3));
-  }
-
-  Ndir.Multiply(myOffset / R);
-  // P(u)
-  theValue.ChangeCoord().Add(Ndir);
-  // P'(u)
-  theD1.Add(gp_Vec2d(DNdir));
-}
-
-void Geom2dEvaluator_OffsetCurve::CalculateD2(      gp_Pnt2d& theValue,
-                                                    gp_Vec2d& theD1,
-                                                    gp_Vec2d& theD2,
-                                              const gp_Vec2d& theD3,
-                                              const Standard_Boolean theIsDirChange) const
-{
-  // P(u) = p(u) + Offset * Ndir / R
-  // with R = || p' ^ Z|| and Ndir = P' ^ Z
-
-  // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
-
-  // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
-  //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
-
-  gp_XY Ndir(theD1.Y(), -theD1.X());
-  gp_XY DNdir(theD2.Y(), -theD2.X());
-  gp_XY D2Ndir(theD3.Y(), -theD3.X());
-  Standard_Real R2 = Ndir.SquareModulus();
-  Standard_Real R = Sqrt(R2);
-  Standard_Real R3 = R2 * R;
-  Standard_Real R4 = R2 * R2;
-  Standard_Real R5 = R3 * R2;
-  Standard_Real Dr = Ndir.Dot(DNdir);
-  Standard_Real D2r = Ndir.Dot(D2Ndir) + DNdir.Dot(DNdir);
-  if (R5 <= gp::Resolution())
-  {
-    if (R4 <= gp::Resolution())
-      throw Standard_NullValue("Geom2dEvaluator_OffsetCurve: Null derivative");
-    //We try another computation but the stability is not very good dixit ISG.
-    // V2 = P" (U) :
-    D2Ndir.Subtract(DNdir.Multiplied(2.0 * Dr / R2));
-    D2Ndir.Add(Ndir.Multiplied(((3.0 * Dr * Dr) / R4) - (D2r / R2)));
-    D2Ndir.Multiply(myOffset / R);
-
-    // V1 = P' (U) :
-    DNdir.Multiply(R);
-    DNdir.Subtract(Ndir.Multiplied(Dr / R));
-    DNdir.Multiply(myOffset / R2);
-  }
-  else
-  {
-    // Same computation as IICURV in EUCLID-IS because the stability is better.
-    // V2 = P" (U) :
-    D2Ndir.Multiply(myOffset / R);
-    D2Ndir.Subtract(DNdir.Multiplied(2.0 * myOffset * Dr / R3));
-    D2Ndir.Add(Ndir.Multiplied(myOffset * (((3.0 * Dr * Dr) / R5) - (D2r / R3))));
-
-    // V1 = P' (U) 
-    DNdir.Multiply(myOffset / R);
-    DNdir.Subtract(Ndir.Multiplied(myOffset * Dr / R3));
-  }
-
-  Ndir.Multiply(myOffset / R);
-  // P(u)
-  theValue.ChangeCoord().Add(Ndir);
-  // P'(u) :
-  theD1.Add(gp_Vec2d(DNdir));
-  // P"(u) :
-  if (theIsDirChange)
-    theD2.Reverse();
-  theD2.Add(gp_Vec2d(D2Ndir));
-}
-
-void Geom2dEvaluator_OffsetCurve::CalculateD3(      gp_Pnt2d& theValue,
-                                                    gp_Vec2d& theD1,
-                                                    gp_Vec2d& theD2,
-                                                    gp_Vec2d& theD3,
-                                              const gp_Vec2d& theD4,
-                                              const Standard_Boolean theIsDirChange) const
-{
-  // P(u) = p(u) + Offset * Ndir / R
-  // with R = || p' ^ Z|| and Ndir = P' ^ Z
-
-  // P'(u)  = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
-
-  // P"(u)  = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
-  //          Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
-
-  // P"'(u) = p"'(u) + (Offset / R) * (D3Ndir - (3.0 * Dr/R**2 ) * D2Ndir -
-  //          (3.0 * D2r / R2) * DNdir) + (3.0 * Dr * Dr / R4) * DNdir -
-  //          (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir +
-  //          (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir
-
-  gp_XY Ndir(theD1.Y(), -theD1.X());
-  gp_XY DNdir(theD2.Y(), -theD2.X());
-  gp_XY D2Ndir(theD3.Y(), -theD3.X());
-  gp_XY D3Ndir(theD4.Y(), -theD4.X());
-  Standard_Real R2 = Ndir.SquareModulus();
-  Standard_Real R = Sqrt(R2);
-  Standard_Real R3 = R2 * R;
-  Standard_Real R4 = R2 * R2;
-  Standard_Real R5 = R3 * R2;
-  Standard_Real R6 = R3 * R3;
-  Standard_Real R7 = R5 * R2;
-  Standard_Real Dr = Ndir.Dot(DNdir);
-  Standard_Real D2r = Ndir.Dot(D2Ndir) + DNdir.Dot(DNdir);
-  Standard_Real D3r = Ndir.Dot(D3Ndir) + 3.0 * DNdir.Dot(D2Ndir);
-
-  if (R7 <= gp::Resolution())
-  {
-    if (R6 <= gp::Resolution())
-      throw Standard_NullValue("Geom2dEvaluator_OffsetCurve: Null derivative");
-    //We try another computation but the stability is not very good dixit ISG.
-    // V3 = P"' (U) :
-    D3Ndir.Subtract(D2Ndir.Multiplied(3.0 * myOffset * Dr / R2));
-    D3Ndir.Subtract(
-      (DNdir.Multiplied((3.0 * myOffset) * ((D2r / R2) + (Dr*Dr) / R4))));
-    D3Ndir.Add(Ndir.Multiplied(
-      (myOffset * (6.0*Dr*Dr / R4 + 6.0*Dr*D2r / R4 - 15.0*Dr*Dr*Dr / R6 - D3r))));
-    D3Ndir.Multiply(myOffset / R);
-    // V2 = P" (U) :
-    R4 = R2 * R2;
-    D2Ndir.Subtract(DNdir.Multiplied(2.0 * Dr / R2));
-    D2Ndir.Subtract(Ndir.Multiplied(((3.0 * Dr * Dr) / R4) - (D2r / R2)));
-    D2Ndir.Multiply(myOffset / R);
-    // V1 = P' (U) :
-    DNdir.Multiply(R);
-    DNdir.Subtract(Ndir.Multiplied(Dr / R));
-    DNdir.Multiply(myOffset / R2);
-  }
-  else
-  {
-    // Same computation as IICURV in EUCLID-IS because the stability is better.
-    // V3 = P"' (U) :
-    D3Ndir.Multiply(myOffset / R);
-    D3Ndir.Subtract(D2Ndir.Multiplied(3.0 * myOffset * Dr / R3));
-    D3Ndir.Subtract(DNdir.Multiplied(
-      ((3.0 * myOffset) * ((D2r / R3) + (Dr*Dr) / R5))));
-    D3Ndir.Add(Ndir.Multiplied(
-      (myOffset * (6.0*Dr*Dr / R5 + 6.0*Dr*D2r / R5 - 15.0*Dr*Dr*Dr / R7 - D3r))));
-    // V2 = P" (U) :
-    D2Ndir.Multiply(myOffset / R);
-    D2Ndir.Subtract(DNdir.Multiplied(2.0 * myOffset * Dr / R3));
-    D2Ndir.Subtract(Ndir.Multiplied(
-      myOffset * (((3.0 * Dr * Dr) / R5) - (D2r / R3))));
-    // V1 = P' (U) :
-    DNdir.Multiply(myOffset / R);
-    DNdir.Subtract(Ndir.Multiplied(myOffset * Dr / R3));
-  }
-
-  Ndir.Multiply(myOffset / R);
-  // P(u)
-  theValue.ChangeCoord().Add(Ndir);
-  // P'(u) :
-  theD1.Add(gp_Vec2d(DNdir));
-  // P"(u)
-  theD2.Add(gp_Vec2d(D2Ndir));
-  // P"'(u)
-  if (theIsDirChange)
-    theD3.Reverse();
-  theD3.Add(gp_Vec2d(D2Ndir));
-}
 
 
 Standard_Boolean Geom2dEvaluator_OffsetCurve::AdjustDerivative(