0023706: Cannot project point on curve
authornbv <nbv@opencascade.com>
Thu, 13 Jun 2013 11:12:06 +0000 (15:12 +0400)
committernbv <nbv@opencascade.com>
Thu, 13 Jun 2013 11:12:06 +0000 (15:12 +0400)
1.   Approximation of derivative (by Taylor-series and by three points).
2.   Some methods (Degree(), GetType(), D0(), D3(), DN()) are added.
3.   Getting of subInterval's boundaries.
4.   Algorithm for checking if 1st derivative is equal to zero is amended.
5.   Cases are controlled when extrema or Project point do not exist.
6.   GetNormal() function for gp_Vec2d was added.
7.   Computing of Value, D0, D1, D2 and D3 for offset curves was changed.
8.   Limitation of tolerance for derivative computing was added.
9.   Methods for computing trihedron in singularity point are added.
10. Test tests/bugs/moddata_3/bug23706 is added.
11. Restriction on the LastParameter for visualization of 3-D curves. Calling PlotCurve(...) function for last interval.
12. LProp package is modified for tangent computing in singularity point (LProp_CLProps, LProp_SLProps).
13. Added test cases for issue.
Deleting bad test cases for this fix

93 files changed:
src/AdvApprox/AdvApprox_ApproxAFunction.cxx
src/Approx/Approx_SweepApproximation.cxx
src/DrawTrSurf/DrawTrSurf_Drawable.cxx
src/Extrema/Extrema_CurveTool.cdl
src/Extrema/Extrema_CurveTool.lxx
src/Extrema/Extrema_FuncExtCC.cdl
src/Extrema/Extrema_FuncExtCC.gxx
src/Extrema/Extrema_FuncExtCC.lxx
src/Extrema/Extrema_FuncExtPC.cdl
src/Extrema/Extrema_FuncExtPC.gxx
src/Extrema/Extrema_GExtCC2d.gxx
src/Extrema/Extrema_GExtPC.gxx
src/Extrema/Extrema_GenExtCC.gxx
src/Extrema/Extrema_GenExtPC.gxx
src/Geom/Geom_OffsetCurve.cxx
src/Geom2d/Geom2d_OffsetCurve.cxx
src/Geom2dInt/Geom2dInt_CurveTool.cdl
src/Geom2dInt/Geom2dInt_CurveTool.lxx
src/GeomFill/GeomFill_Frenet.cdl
src/GeomFill/GeomFill_Frenet.cxx
src/GeomFill/GeomFill_SnglrFunc.cdl
src/GeomFill/GeomFill_SnglrFunc.cxx
src/GeometryTest/GeometryTest_APICommands.cxx
src/GeometryTest/GeometryTest_CurveCommands.cxx
src/GeomliteTest/GeomliteTest_API2dCommands.cxx
src/HLRBRep/HLRBRep_CurveTool.cdl
src/HLRBRep/HLRBRep_CurveTool.lxx
src/IntCurve/IntCurve_IntCurveCurveGen.gxx
src/IntCurve/IntCurve_UserIntConicCurveGen.gxx
src/LProp/LProp_CLProps.cdl
src/LProp/LProp_CLProps.gxx
src/LProp/LProp_SLProps.cdl
src/LProp/LProp_SLProps.gxx
src/gp/gp_Vec.cdl
src/gp/gp_Vec2d.cdl
src/gp/gp_Vec2d.lxx
tests/bugs/modalg_5/bug23706_1 [new file with mode: 0755]
tests/bugs/modalg_5/bug23706_10 [new file with mode: 0755]
tests/bugs/modalg_5/bug23706_11 [new file with mode: 0755]
tests/bugs/modalg_5/bug23706_12 [new file with mode: 0755]
tests/bugs/modalg_5/bug23706_13 [new file with mode: 0755]
tests/bugs/modalg_5/bug23706_14 [new file with mode: 0755]
tests/bugs/modalg_5/bug23706_15 [new file with mode: 0755]
tests/bugs/modalg_5/bug23706_16 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_17 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_19 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_2 [new file with mode: 0755]
tests/bugs/modalg_5/bug23706_20 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_21 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_22 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_24 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_26 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_27 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_28 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_29 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_3 [new file with mode: 0755]
tests/bugs/modalg_5/bug23706_31 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_32 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_33 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_34 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_36 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_37 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_38 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_39 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_4 [new file with mode: 0755]
tests/bugs/modalg_5/bug23706_40 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_41 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_42 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_43 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_44 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_45 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_46 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_47 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_48 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_49 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_5 [new file with mode: 0755]
tests/bugs/modalg_5/bug23706_50 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_51 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_52 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_53 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_54 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_55 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_56 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_57 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_58 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_59 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_6 [new file with mode: 0755]
tests/bugs/modalg_5/bug23706_60 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_61 [new file with mode: 0644]
tests/bugs/modalg_5/bug23706_7 [new file with mode: 0755]
tests/bugs/modalg_5/bug23706_8 [new file with mode: 0755]
tests/bugs/modalg_5/bug23706_9 [new file with mode: 0755]
tests/bugs/moddata_3/bug23706 [new file with mode: 0644]

index 5bbecf7..af57b2a 100755 (executable)
@@ -786,20 +786,20 @@ void AdvApprox_ApproxAFunction::Perform(const Standard_Integer Num1DSS,
 
   TColStd_Array1OfReal AverageError(1,myMaxSegments * TotalNumSS) ;
   
-  Approximation (TotalDimension,  TotalNumSS, LocalDimension,  
-                myFirst,         myLast,
-                *(AdvApprox_EvaluatorFunction*)myEvaluator, 
-                 CutTool,    
-                ContinuityOrder,
-                NumMaxCoeffs,    myMaxSegments,
-                LocalTolerances, code_precis,
-                 NumCurves,                           // Nombre de courbe en sortie
-                NumCoeffPerCurvePtr->ChangeArray1(), // Nbre de coeff par courbe
-                LocalCoefficientsPtr->ChangeArray1(),// Les Coeffs solutions
-                IntervalsPtr->ChangeArray1(),        // La Table des decoupes
-                ErrorMax,                            // Majoration de l'erreur
-                AverageError,                        // Erreur moyenne constatee
-                ErrorCode) ;
+  Approximation ( TotalDimension,
+                  TotalNumSS,   LocalDimension,
+                  myFirst,      myLast,
+                  *(AdvApprox_EvaluatorFunction*)myEvaluator, 
+                  CutTool, ContinuityOrder, NumMaxCoeffs,
+                  myMaxSegments, LocalTolerances, code_precis,
+                  NumCurves,                            // Nombre de courbe en sortie
+                  NumCoeffPerCurvePtr->ChangeArray1(),  // Nbre de coeff par courbe
+                  LocalCoefficientsPtr->ChangeArray1(), // Les Coeffs solutions
+                  IntervalsPtr->ChangeArray1(),         // La Table des decoupes
+                  ErrorMax,                             // Majoration de l'erreur
+                  AverageError,                         // Erreur moyenne constatee
+                  ErrorCode) ;
+  
   if (ErrorCode == 0 || ErrorCode == -1)    {
     //
     // si tout est OK ou bien on a un resultat dont l une des erreurs max est
index ff478da..2743fd4 100755 (executable)
@@ -478,11 +478,11 @@ Standard_Boolean Approx_SweepApproximation::D1(const Standard_Real Param,
       myPoles->ChangeValue(ii).ChangeCoord()
        -= Translation.XYZ();
       // Homothety on all. 
-      myDPoles->ChangeValue(ii) *= myWeigths->Value(ii);
+      const Standard_Real aWeight = myWeigths->Value(ii);
+      myDPoles->ChangeValue(ii) *= aWeight;
       Vaux.SetXYZ( myPoles->Value(ii).Coord());
       myDPoles->ChangeValue(ii) += myDWeigths->Value(ii)*Vaux;
-      myPoles->ChangeValue(ii).ChangeCoord() 
-       *= myWeigths->Value(ii); // for the cash
+      myPoles->ChangeValue(ii).ChangeCoord() *= aWeight; // for the cash
     }
       
 
index 672bd43..3a966e7 100755 (executable)
@@ -127,21 +127,25 @@ static void PlotCurve (Draw_Display& aDisplay,
 //=======================================================================
 
 void DrawTrSurf_Drawable::DrawCurveOn (Adaptor3d_Curve&   C, 
-                                      Draw_Display& aDisplay) const
+                                       Draw_Display& aDisplay) const
 {
   gp_Pnt P;
-  if (myDrawMode == 1) {
+  if (myDrawMode == 1)
+  {
     Standard_Real Fleche = myDeflection/aDisplay.Zoom();
     GCPnts_UniformDeflection LineVu(C,Fleche);
-    if (LineVu.IsDone()) {
+    if (LineVu.IsDone())
+    {
       aDisplay.MoveTo(LineVu.Value(1));
-      for (Standard_Integer i = 2; i <= LineVu.NbPoints(); i++) {
-       aDisplay.DrawTo(LineVu.Value(i));
+      for (Standard_Integer i = 2; i <= LineVu.NbPoints(); i++)
+      {
+        aDisplay.DrawTo(LineVu.Value(i));
       }
-    }  
+    }
   }
-  else {
-    Standard_Real j;
+  else
+  {
+    Standard_Integer j;
     Standard_Integer intrv, nbintv = C.NbIntervals(GeomAbs_CN);
     TColStd_Array1OfReal TI(1,nbintv+1);
     C.Intervals(TI,GeomAbs_CN);
@@ -150,36 +154,44 @@ void DrawTrSurf_Drawable::DrawCurveOn (Adaptor3d_Curve&   C,
     GeomAbs_CurveType CurvType = C.GetType();
     gp_Pnt aPPnt=P, aNPnt;
 
-    for (intrv = 1; intrv <= nbintv; intrv++) {
+    for (intrv = 1; intrv <= nbintv; intrv++)
+    {
       Standard_Real t = TI(intrv);
       Standard_Real step = (TI(intrv+1) - t) / myDiscret;
 
-      switch (CurvType) {
-      case GeomAbs_Line :
-       break;
-      case GeomAbs_Circle :
-      case GeomAbs_Ellipse :
-       for (j = 1; j < myDiscret; j++) {
-         t += step;
-         C.D0(t,P);
-         aDisplay.DrawTo(P);
-       }
-       break;
-      case GeomAbs_Parabola :
-      case GeomAbs_Hyperbola :
-      case GeomAbs_BezierCurve :
-      case GeomAbs_BSplineCurve :
-      case GeomAbs_OtherCurve :
-       for (j = 1; j <= myDiscret/2; j++) {
-         C.D0 (t+step*2., aNPnt);
+      switch (CurvType)
+      {
+      case GeomAbs_Line:
+        break;
+      case GeomAbs_Circle:
+      case GeomAbs_Ellipse:
+        for (j = 1; j < myDiscret; j++)
+        {
+          t += step;
+          C.D0(t,P);
+          aDisplay.DrawTo(P);
+        }
+        break;
+      case GeomAbs_Parabola:
+      case GeomAbs_Hyperbola:
+      case GeomAbs_BezierCurve:
+      case GeomAbs_BSplineCurve:
+      case GeomAbs_OtherCurve:
+        const Standard_Integer nIter = myDiscret/2;
+        for (j = 1; j < nIter; j++)
+        {
+          const Standard_Real t1 = t+step*2.;
+          C.D0 (t1, aNPnt);
           PlotCurve (aDisplay, C, t, step, aPPnt, aNPnt);
-         aPPnt = aNPnt;
-         t += step*2.;
-       }
-       break;
+          aPPnt = aNPnt;
+          t = t1;
+        }
+
+        break;
       }
 
       C.D0(TI(intrv+1),P);
+      PlotCurve (aDisplay, C, t, step, aPPnt, P);
       aDisplay.DrawTo(P);
     }
   }
index 12c5f66..13e6726 100755 (executable)
@@ -83,6 +83,11 @@ is
                    U : Real    from Standard)
     returns Pnt from gp;
        ---C++: inline
+
+    D0 (myclass; C :     Curve from Adaptor3d;
+                 U :     Real    from Standard;
+                 P : out Pnt from gp);
+       ---C++: inline
     
     D1 (myclass; C :     Curve from Adaptor3d;
                  U :     Real    from Standard;
@@ -90,10 +95,22 @@ is
                  V : out Vec from gp);
        ---C++: inline
     
-    D2 (myclass; C      :     Curve from Adaptor3d;
+    D2 (myclass; C          : Curve from Adaptor3d;
+                 U          : Real  from Standard;
+                 P          : out Pnt from gp;
+                 V1, V2     : out Vec from gp);
+       ---C++: inline
+
+    D3 (myclass; C          : Curve from Adaptor3d;
+                 U          : Real  from Standard;
+                 P          : out Pnt from gp;
+                 V1, V2, V3 : out Vec from gp);
+       ---C++: inline
+
+    DN (myclass; C      :     Curve from Adaptor3d;
                  U      :     Real    from Standard;
-                 P      : out Pnt from gp;
-                 V1, V2 : out Vec from gp);
+                 N      :     Integer from Standard)
+       returns Vec from gp;          
        ---C++: inline
 
    Line(myclass; C : Curve from Adaptor3d) returns Lin from gp;
index 6f84904..85a6299 100755 (executable)
@@ -62,6 +62,17 @@ inline gp_Pnt Extrema_CurveTool::Value(const Adaptor3d_Curve& C,
   return C.Value(U);
 }
 
+//=======================================================================
+//function : D0
+//purpose  : 
+//=======================================================================
+
+inline void Extrema_CurveTool::D0(const Adaptor3d_Curve& C, 
+                                       const Standard_Real U, 
+                                       gp_Pnt& P)
+{
+  C.D0(U, P);
+}
 
 //=======================================================================
 //function : D1
@@ -90,6 +101,31 @@ inline void Extrema_CurveTool::D2(const Adaptor3d_Curve& C,
   C.D2(U, P, V1, V2);
 }
 
+//=======================================================================
+//function : D3
+//purpose  : 
+//=======================================================================
+
+ inline void Extrema_CurveTool::D3(const Adaptor3d_Curve& C, 
+                                        const Standard_Real U, 
+                                        gp_Pnt& P, 
+                                        gp_Vec& V1, 
+                                        gp_Vec& V2, 
+                                        gp_Vec& V3)
+{
+  C.D3(U, P, V1, V2, V3);
+}
+
+//=======================================================================
+//function : DN
+//purpose  : 
+//=======================================================================
+inline gp_Vec Extrema_CurveTool::DN(const Adaptor3d_Curve& C, 
+                                            const Standard_Real U, 
+                                            const Standard_Integer N)
+{
+  return C.DN(U, N);
+}
 
 
 //=======================================================================
index e49d6b7..45ebdcc 100755 (executable)
@@ -50,7 +50,6 @@ is
     ---Purpose:
 
     SetCurve (me: in out; theRank: Integer; C1: Curve1);
-    ---C++: inline
     ---Purpose:
 
     SetTolerance (me: in out; theTol: Real);
@@ -103,6 +102,14 @@ is
         ---Purpose: Returns a tolerance specified in the constructor
         --          or in SetTolerance() method.
 
+    SubIntervalInitialize(me: in out; theUfirst, theUlast: Vector);
+       ---Purpose: Determines of boundaries of subinterval for find of root.
+    
+    SearchOfTolerance(me: in out; C: Address from Standard) returns Real from Standard;
+       ---Purpose: Computes a Tol value. If 1st derivative of curve
+       --          |D1|<Tol, it is considered D1=0.      
+
+
 fields
     myC1    : Address from Standard;
     myC2    : Address from Standard;
@@ -117,4 +124,14 @@ fields
     mySqDist: SequenceOfReal from TColStd;
     myPoints: SeqPOnC;
 
+    myTolC1,myTolC2: Real from Standard;  -- toolerance for derivate
+
+--Supremum of search 1st non-zero derivative    
+    myMaxDerivOrderC1, myMaxDerivOrderC2: Integer from Standard;
+    
+--boundaries of subinterval for find of root
+    myUinfium, myUsupremum: Real from Standard; -- C1 curve
+    myVinfium, myVsupremum: Real from Standard; -- C2 curve
+
+
 end FuncExtCC;
index 8fae8b8..7eda01c 100755 (executable)
@@ -38,126 +38,690 @@ et F2 sont egales a:
   { Dvf2(u,v) = ||Dv|| + C2C1.Dvv/||Dv||- F2(u,v)*Dv*Dvv/||Dv||**2
 
 ----------------------------------------------------------------------------*/
-//=============================================================================
 
-#define Tol  1.e-20
-#define delta 1.e-9
+#include <Precision.hxx>
+
+
+static const Standard_Real MinTol    = 1.e-20;
+static const Standard_Real TolFactor = 1.e-12;
+static const Standard_Real MinStep   = 1e-7;
+
+static const Standard_Integer MaxOrder = 3;
+
+
+
+//=============================================================================
+Standard_Real Extrema_FuncExtCC::SearchOfTolerance(const Standard_Address C)
+  {
+  const Standard_Integer NPoint = 10;
+  Standard_Real aStartParam, anEndParam;
+  
+  if(C==myC1)
+    {
+    aStartParam = myUinfium;
+    anEndParam = myUsupremum;
+    }
+  else if(C==myC2)
+    {
+    aStartParam = myVinfium;
+    anEndParam = myVsupremum;
+    }
+  else
+    {
+#ifdef DEB
+    cout << "+++ Function Extrema_FuncExtCC::SearchOfTolerance(...)" << endl;
+    cout << "Warning: No curve for tolerance computing!---"<<endl;
+#endif
+    return MinTol;
+    }
+    
+  const Standard_Real aStep = (anEndParam - aStartParam)/(Standard_Real)NPoint;
+    
+  Standard_Integer aNum = 0;
+  Standard_Real aMax = -Precision::Infinite();  //Maximum value of 1st derivative
+                                                //(it is computed with using NPoint point)
+    
+  do
+    {
+    Standard_Real u = aStartParam + aNum*aStep;    //parameter for every point
+    if(u > anEndParam)
+      u = anEndParam;
+    
+    Pnt Ptemp;  //empty point (is not used below)
+    Vec VDer;   // 1st derivative vector
+    Tool1::D1(*((Curve1*)C), u, Ptemp, VDer);
+    Standard_Real vm = VDer.Magnitude();
+    if(vm > aMax)
+      aMax = vm;
+    }
+  while(++aNum < NPoint+1);
+  
+  return Max(aMax*TolFactor,MinTol);
+  }
 
+//=============================================================================
 Extrema_FuncExtCC::Extrema_FuncExtCC(const Standard_Real thetol) : myC1 (0), myC2 (0), myTol (thetol)
-{
-}
+  {
+  math_Vector V1(1,2), V2(1,2);
+  V1(1) = 0.0;
+  V2(1) = 0.0;
+  V1(2) = 0.0;
+  V2(2) = 0.0;
+  SubIntervalInitialize(V1, V2);
+  myMaxDerivOrderC1 = 0;
+  myTolC1=MinTol;
+  myMaxDerivOrderC2 = 0;
+  myTolC2=MinTol;
+  }
 
+//=============================================================================
 Extrema_FuncExtCC::Extrema_FuncExtCC (const Curve1& C1,
                                      const Curve2& C2,
                                      const Standard_Real thetol) :
                                        myC1 ((Standard_Address)&C1), myC2 ((Standard_Address)&C2),
                                        myTol (thetol)
-{
-}
+  {
+  math_Vector V1(1,2), V2(1,2);
+  V1(1) = Tool1::FirstParameter(*((Curve1*)myC1));
+  V2(1) = Tool1::LastParameter(*((Curve1*)myC1));
+  V1(2) = Tool2::FirstParameter(*((Curve2*)myC2));
+  V2(2) = Tool2::LastParameter(*((Curve2*)myC2));
+  SubIntervalInitialize(V1, V2);
+                        
+  switch(Tool1::GetType(*((Curve1*)myC1)))
+    {
+    case GeomAbs_BezierCurve:
+    case GeomAbs_BSplineCurve:
+    case GeomAbs_OtherCurve:
+      myMaxDerivOrderC1 = MaxOrder;
+      myTolC1 = SearchOfTolerance((Standard_Address)&C1);
+      break;
+    default:
+      myMaxDerivOrderC1 = 0;
+      myTolC1=MinTol;
+      break;
+    }
+    
+  switch(Tool2::GetType(*((Curve2*)myC2)))
+    {
+    case GeomAbs_BezierCurve:
+    case GeomAbs_BSplineCurve:
+    case GeomAbs_OtherCurve:
+      myMaxDerivOrderC2 = MaxOrder;
+      myTolC2 = SearchOfTolerance((Standard_Address)&C2);
+      break;
+    default:
+      myMaxDerivOrderC2 = 0;
+      myTolC2=MinTol;
+      break;
+    }
+  }
 
-Standard_Boolean Extrema_FuncExtCC::Value (const math_Vector& UV, 
-                                          math_Vector&       F)
-{
+//=============================================================================
+void Extrema_FuncExtCC::SetCurve (const Standard_Integer theRank, const Curve1& C)
+  {
+  Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_FuncExtCC::SetCurve()")
   
+  if (theRank == 1)
+    {
+    myC1 = (Standard_Address)&C;
+    switch(/*Tool1::GetType(*((Curve1*)myC1))*/ C.GetType())
+      {
+      case GeomAbs_BezierCurve:
+      case GeomAbs_BSplineCurve:
+      case GeomAbs_OtherCurve:
+        myMaxDerivOrderC1 = MaxOrder;
+        myTolC1 = SearchOfTolerance((Standard_Address)&C);
+        break;
+      default:
+        myMaxDerivOrderC1 = 0;
+        myTolC1=MinTol;
+        break;
+      }
+    }
+  else if (theRank == 2)
+    {
+    myC2 = (Standard_Address)&C;
+    switch(/*Tool2::GetType(*((Curve2*)myC2))*/C.GetType())
+      {
+      case GeomAbs_BezierCurve:
+      case GeomAbs_BSplineCurve:
+      case GeomAbs_OtherCurve:
+        myMaxDerivOrderC2 = MaxOrder;
+        myTolC2 = SearchOfTolerance((Standard_Address)&C);
+        break;
+      default:
+        myMaxDerivOrderC2 = 0;
+        myTolC2=MinTol;
+        break;
+      }
+    }
+  }
+
+//=============================================================================
+Standard_Boolean Extrema_FuncExtCC::Value (const math_Vector& UV, math_Vector& F)
+  {
   myU = UV(1);
   myV = UV(2);
   Tool1::D1(*((Curve1*)myC1), myU,myP1,myDu);
   Tool2::D1(*((Curve2*)myC2), myV,myP2,myDv);
+
   Vec P1P2 (myP1,myP2);
 
   Standard_Real Ndu = myDu.Magnitude();
-  if (Ndu <= Tol) {
-    Pnt P1, P2;
-    P1 = Tool1::Value(*((Curve1*)myC1), myU-delta);
-    P2 = Tool1::Value(*((Curve1*)myC1), myU+delta);
-    Vec V(P1,P2);
-    myDu = V;
-    Ndu = myDu.Magnitude();
-    if (Ndu <= Tol) {
-      return Standard_False;
+
+
+  if(myMaxDerivOrderC1 != 0)
+    {
+    if (Ndu <= myTolC1)
+      {
+    //Derivative is approximated by Taylor-series
+      const Standard_Real DivisionFactor = 1.e-3;
+      Standard_Real du;
+      if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
+        du = 0.0;
+      else
+        du = myUsupremum-myUinfium;
+        
+      const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
+      
+      Standard_Integer n = 1; //Derivative order
+      Vec V;
+      Standard_Boolean IsDeriveFound;
+      
+      do
+        {
+        V = Tool1::DN(*((Curve1*)myC1),myU,++n);
+        Ndu = V.Magnitude();
+        IsDeriveFound = (Ndu > myTolC1);
+        }
+      while(!IsDeriveFound && n < myMaxDerivOrderC1);
+      
+      if(IsDeriveFound)
+        {
+        Standard_Real u;
+        
+        if(myU-myUinfium < aDelta)
+          u = myU+aDelta;
+        else
+          u = myU-aDelta;
+          
+        Pnt P1, P2;
+        Tool1::D0(*((Curve1*)myC1),Min(myU, u),P1);
+        Tool1::D0(*((Curve1*)myC1),Max(myU, u),P2);
+        
+        Vec V1(P1,P2);
+        Standard_Real aDirFactor = V.Dot(V1);
+        
+        if(aDirFactor < 0.0)
+          myDu = -V;
+        else
+          myDu = V;
+        }//if(IsDeriveFound)
+      else
+        {
+  //Derivative is approximated by three points
+
+        Pnt Ptemp; //(0,0,0)-coordinate
+        Pnt P1, P2, P3;
+        Standard_Boolean IsParameterGrown;
+                  
+        if(myU-myUinfium < 2*aDelta)
+          {
+          Tool1::D0(*((Curve1*)myC1),myU,P1);
+          Tool1::D0(*((Curve1*)myC1),myU+aDelta,P2);
+          Tool1::D0(*((Curve1*)myC1),myU+2*aDelta,P3);
+          IsParameterGrown = Standard_True;
+          }
+        else
+          {
+          Tool1::D0(*((Curve1*)myC1),myU-2*aDelta,P1);
+          Tool1::D0(*((Curve1*)myC1),myU-aDelta,P2);
+          Tool1::D0(*((Curve1*)myC1),myU,P3);
+          IsParameterGrown = Standard_False;
+          }
+          
+        Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
+        
+        if(IsParameterGrown)
+          myDu=-3*V1+4*V2-V3;
+        else
+          myDu=V1-4*V2+3*V3;
+        }//else of if(IsDeriveFound)
+      Ndu = myDu.Magnitude();      
+      }//if (Ndu <= myTolC1) condition
+    }//if(myMaxDerivOrder != 0)
+
+  if (Ndu <= MinTol)
+    {
+#ifdef DEB
+    cout << "+++Function Extrema_FuncExtCC::Value(...)." << endl;
+    cout << "Warning: 1st derivative of C1 is equal to zero!---"<<endl;
+#endif
+    return Standard_False;
     }
-  }
+
 
   Standard_Real Ndv = myDv.Magnitude();
-  if (Ndv <= Tol) { 
-    // Traitement des singularite, on approche la Tangente
-    // par une corde
-    Pnt P1, P2;
-    P1 = Tool2::Value(*((Curve2*)myC2), myV-delta);
-    P2 = Tool2::Value(*((Curve2*)myC2), myV+delta);
-    Vec V(P1,P2);
-    myDv = V;
-    Ndv = myDv.Magnitude();
-    if (Ndv <= Tol) {
-      return Standard_False;
+
+  if(myMaxDerivOrderC2 != 0)
+    {
+    if (Ndv <= myTolC2)
+      { 
+      const Standard_Real DivisionFactor = 1.e-3;
+      Standard_Real dv;
+      if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
+        dv = 0.0;
+      else
+        dv = myVsupremum-myVinfium;
+        
+      const Standard_Real aDelta = Max(dv*DivisionFactor,MinStep);
+
+//Derivative is approximated by Taylor-series
+      
+      Standard_Integer n = 1; //Derivative order
+      Vec V;
+      Standard_Boolean IsDeriveFound;
+      
+      do
+        {
+        V = Tool2::DN(*((Curve2*)myC2),myV,++n);
+        Ndv = V.Magnitude();
+        IsDeriveFound = (Ndv > myTolC2);
+        }
+      while(!IsDeriveFound && n < myMaxDerivOrderC2);
+      
+      if(IsDeriveFound)
+        {
+        Standard_Real v;
+        
+        if(myV-myVinfium < aDelta)
+          v = myV+aDelta;
+        else
+          v = myV-aDelta;
+          
+        Pnt P1, P2;
+        Tool2::D0(*((Curve2*)myC2),Min(myV, v),P1);
+        Tool2::D0(*((Curve2*)myC2),Max(myV, v),P2);
+        
+        Vec V1(P1,P2);
+        Standard_Real aDirFactor = V.Dot(V1);
+        
+        if(aDirFactor < 0.0)
+          myDv = -V;
+        else
+          myDv = V;
+        }//if(IsDeriveFound)
+      else
+        {
+  //Derivative is approximated by three points
+
+        Pnt Ptemp; //(0,0,0)-coordinate
+        Pnt P1, P2, P3;
+        Standard_Boolean IsParameterGrown;
+                  
+        if(myV-myVinfium < 2*aDelta)
+          {
+          Tool2::D0(*((Curve2*)myC2),myV,P1);
+          Tool2::D0(*((Curve2*)myC2),myV+aDelta,P2);
+          Tool2::D0(*((Curve2*)myC2),myV+2*aDelta,P3);
+          IsParameterGrown = Standard_True;
+          }
+        else
+          {
+          Tool2::D0(*((Curve2*)myC2),myV-2*aDelta,P1);
+          Tool2::D0(*((Curve2*)myC2),myV-aDelta,P2);
+          Tool2::D0(*((Curve2*)myC2),myV,P3);
+          IsParameterGrown = Standard_False;
+          }
+          
+        Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
+        
+        if(IsParameterGrown)
+          myDv=-3*V1+4*V2-V3;
+        else
+          myDv=V1-4*V2+3*V3;
+        }//else of if(IsDeriveFound)
+        
+        Ndv = myDv.Magnitude();
+      }//if (Ndv <= myTolC2)
+    }//if(myMaxDerivOrder != 0)
+
+  if (Ndv <= MinTol)
+    {
+#ifdef DEB
+    cout << "+++Function Extrema_FuncExtCC::Value(...)." << endl;
+    cout << "1st derivative of C2 is equal to zero!---"<<endl;
+#endif
+    return Standard_False;
     }
-  }
 
   F(1) = P1P2.Dot(myDu)/Ndu;
   F(2) = P1P2.Dot(myDv)/Ndv;
   return Standard_True;
-}
+  }
 //=============================================================================
 
 Standard_Boolean Extrema_FuncExtCC::Derivatives (const math_Vector& UV, 
                                                 math_Matrix& Df)
-{
+  {
   math_Vector F(1,2);
   return Values(UV,F,Df);
-}
+  }
 //=============================================================================
+
 Standard_Boolean Extrema_FuncExtCC::Values (const math_Vector& UV, 
-                                           math_Vector& F,
-                                           math_Matrix& Df)
-{
+                                            math_Vector& F,
+                                            math_Matrix& Df)
+  {
   myU = UV(1);
   myV = UV(2);
-  Vec Duu, Dvv;
-  Tool1::D2(*((Curve1*)myC1), myU,myP1,myDu,Duu);
-  Tool2::D2(*((Curve2*)myC2), myV,myP2,myDv,Dvv);
+
+  if(Value(UV, F) == Standard_False) //Computes F, myDu, myDv
+    {
+#ifdef DEB
+    cout << "+++Standard_Boolean Extrema_FuncExtCC::Values(...)." << endl;
+    cout << "Warning: No function value found!---"<<endl;
+#endif
+    return Standard_False;
+    }//if(Value(UV, F) == Standard_False)
+
+  Vec Du, Dv, Duu, Dvv;
+  Tool1::D2(*((Curve1*)myC1), myU,myP1,Du,Duu);
+  Tool2::D2(*((Curve2*)myC2), myV,myP2,Dv,Dvv);
+  
+//Calling of "Value(...)" function change class member values. 
+//After running it is necessary to return to previous values.  
+  const Standard_Real myU_old  = myU,    myV_old = myV;
+  const Pnt           myP1_old = myP1,  myP2_old = myP2;
+  const Vec           myDu_old = myDu,  myDv_old = myDv;
+  
+
+//Attention:  aDelta value must be greater than same value for "Value(...)"
+//            function to avoid of points' collisions.
+
+  const Standard_Real DivisionFactor = 0.01;
+  
+  Standard_Real du;
+  if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
+    du = 0.0;
+  else
+    du = myUsupremum-myUinfium;
+    
+  const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep);
+  
+  Standard_Real dv;
+  if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
+    dv = 0.0;
+  else
+    dv = myVsupremum-myVinfium;
+    
+  const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep); 
 
   Vec P1P2 (myP1,myP2);
+  
+  if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
+    {
+//Derivative is approximated by three points
 
-  Standard_Real Ndu = myDu.Magnitude();
-  if (Ndu <= Tol) {
-      Pnt P1, P2;
-    Vec V1;
-    Tool1::D1(*((Curve1*)myC1),myU+delta, P2, Duu);
-    Tool1::D1(*((Curve1*)myC1),myU-delta, P1, V1);
-    Vec V(P1,P2);
-    myDu = V;
-    Duu -= V1;
-    Ndu = myDu.Magnitude();
-    if (Ndu <= Tol) {
-      return Standard_False;
-    }  
-  }
+    math_Vector FF1(1,2), FF2(1,2), FF3(1,2);
+    Standard_Real F1, F2, F3;
 
-  Standard_Real Ndv = myDv.Magnitude();
- if (Ndv <= Tol) {
-    Pnt P1, P2;
-    Vec V1;
-    Tool2::D1(*((Curve2*)myC2),myV+delta, P2, Dvv);
-    Tool2::D1(*((Curve2*)myC2),myV-delta, P1, V1);
-    Vec V(P1,P2);
-    myDv = V;
-    Dvv -= V1;
-    Ndv = myDv.Magnitude();
-    if (Ndv <= Tol) {
-      return Standard_False;
-    }  
-  }
+/////////////////////////// Search of DF1_u derivative (begin) ///////////////////
+    if(myU-myUinfium < 2*aDeltaU)
+      {
+      F1=F(1);
+      math_Vector UV2(1,2), UV3(1,2);
+      UV2(1)=myU+aDeltaU;
+      UV2(2)=myV;
+      UV3(1)=myU+2*aDeltaU;
+      UV3(2)=myV;
+      if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
+        {
+#ifdef DEB
+        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
+        cout << "There are many points close to singularity points "
+            "and which have zero-derivative." << endl;
+        cout << "Try to decrease aDelta variable's value. ---" << endl;
+#endif
+        return Standard_False;
+        }
+        
+      F2 = FF2(1);
+      F3 = FF3(1);
+      
+      Df(1,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU);
+      }//if(myU-myUinfium < 2*aDeltaU)
+    else
+      {
+      F3 = F(1);
+      math_Vector UV2(1,2), UV1(1,2);
+      UV2(1)=myU-aDeltaU;
+      UV2(2)=myV;
+      UV1(1)=myU-2*aDeltaU;
+      UV1(2)=myV;
+
+      if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
+        {
+#ifdef DEB
+        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
+        cout << "There are many points close to singularity points "
+            "and which have zero-derivative." << endl;
+        cout << "Try to decrease aDelta variable's value. ---" << endl;
+#endif
+        return Standard_False;
+        }
+        
+      F1 = FF1(1);
+      F2 = FF2(1);
+
+      Df(1,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU);
+      }//else of if(myU-myUinfium < 2*aDeltaU) condition
+/////////////////////////// Search of DF1_u derivative (end) ///////////////////
+      
+      //Return to previous values
+      myU = myU_old;
+      myV = myV_old;
+
+/////////////////////////// Search of DF1_v derivative (begin) ///////////////////
+    if(myV-myVinfium < 2*aDeltaV)
+      {
+      F1=F(1);
+      math_Vector UV2(1,2), UV3(1,2);
+      UV2(1)=myU;
+      UV2(2)=myV+aDeltaV;
+      UV3(1)=myU;
+      UV3(2)=myV+2*aDeltaV;
+
+      if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
+        {
+#ifdef DEB
+        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
+        cout << "There are many points close to singularity points "
+            "and which have zero-derivative." << endl;
+        cout << "Try to decrease aDelta variable's value. ---" << endl;
+#endif
+        return Standard_False;
+        }
+      F2 = FF2(1);
+      F3 = FF3(1);
+      
+      Df(1,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV);
+      }//if(myV-myVinfium < 2*aDeltaV)
+    else
+      {
+      F3 = F(1);
+      math_Vector UV2(1,2), UV1(1,2);
+      UV2(1)=myU;
+      UV2(2)=myV-aDeltaV;
+      UV1(1)=myU;
+      UV1(2)=myV-2*aDeltaV;
+      if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
+        {
+#ifdef DEB
+        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
+        cout << "There are many points close to singularity points "
+            "and which have zero-derivative." << endl;
+        cout << "Try to decrease aDelta variable's value. ---" << endl;
+#endif
+        return Standard_False;
+        }
+        
+      F1 = FF1(1);
+      F2 = FF2(1);
+
+      Df(1,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV);
+      }//else of if(myV-myVinfium < 2*aDeltaV)
+/////////////////////////// Search of DF1_v derivative (end) ///////////////////
+    
+    //Return to previous values
+    myU = myU_old;
+    myV = myV_old;
+    myP1 = myP1_old,  myP2 = myP2_old;
+    myDu = myDu_old,  myDv = myDv_old;
+    }//if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
+  else
+    {
+    const Standard_Real Ndu = myDu.Magnitude();
+    Df(1,1) = - Ndu + (P1P2.Dot(Duu)/Ndu) - F(1)*(myDu.Dot(Duu)/(Ndu*Ndu));
+    Df(1,2) = myDv.Dot(myDu)/Ndu;
+    }//else of if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
   
-  F(1) = P1P2.Dot(myDu)/Ndu;
-  F(2) = P1P2.Dot(myDv)/Ndv;
+  if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
+    {
+//Derivative is approximated by three points
+
+    math_Vector FF1(1,2), FF2(1,2), FF3(1,2);
+    Standard_Real F1, F2, F3;
+
+/////////////////////////// Search of DF2_v derivative (begin) ///////////////////
+    if(myV-myVinfium < 2*aDeltaV)
+      {
+      F1=F(2);
+      math_Vector UV2(1,2), UV3(1,2);
+      UV2(1)=myU;
+      UV2(2)=myV+aDeltaV;
+      UV3(1)=myU;
+      UV3(2)=myV+2*aDeltaV;
+
+      if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
+        {
+#ifdef DEB
+        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
+        cout << "There are many points close to singularity points "
+            "and which have zero-derivative." << endl;
+        cout << "Try to decrease aDelta variable's value. ---" << endl;
+#endif
+        return Standard_False;
+        }
+        
+      F2 = FF2(2);
+      F3 = FF3(2);
+      
+      Df(2,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV);
+      
+      }//if(myV-myVinfium < 2*aDeltaV)
+    else
+      {
+      F3 = F(2);
+      math_Vector UV2(1,2), UV1(1,2);
+      UV2(1)=myU;
+      UV2(2)=myV-aDeltaV;
+      UV1(1)=myU;
+      UV1(2)=myV-2*aDeltaV;
+
+      if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
+        {
+#ifdef DEB
+        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
+        cout << "There are many points close to singularity points "
+            "and which have zero-derivative." << endl;
+        cout << "Try to decrease aDelta variable's value. ---" << endl;
+#endif
+        return Standard_False;
+        }
+        
+      F1 = FF1(2);
+      F2 = FF2(2);
+
+      Df(2,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV);
+      }//else of if(myV-myVinfium < 2*aDeltaV)
+/////////////////////////// Search of DF2_v derivative (end) ///////////////////
+
+      //Return to previous values
+      myU = myU_old;
+      myV = myV_old;
+
+/////////////////////////// Search of DF2_u derivative (begin) ///////////////////
+    if(myU-myUinfium < 2*aDeltaU)
+      {
+      F1=F(2);
+      math_Vector UV2(1,2), UV3(1,2);
+      UV2(1)=myU+aDeltaU;
+      UV2(2)=myV;
+      UV3(1)=myU+2*aDeltaU;
+      UV3(2)=myV;
+      if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
+        {
+#ifdef DEB
+        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
+        cout << "There are many points close to singularity points "
+            "and which have zero-derivative." << endl;
+        cout << "Try to decrease aDelta variable's value. ---" << endl;
+#endif
+        return Standard_False;
+        }
+        
+      F2 = FF2(2);
+      F3 = FF3(2);
+      
+      Df(2,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU);
+      }//if(myU-myUinfium < 2*aDelta)
+    else
+      {
+      F3 = F(2);
+      math_Vector UV2(1,2), UV1(1,2);
+      UV2(1)=myU-aDeltaU;
+      UV2(2)=myV;
+      UV1(1)=myU-2*aDeltaU;
+      UV1(2)=myV;
+            
+      if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
+        {
+#ifdef DEB
+        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
+        cout << "There are many points close to singularity points "
+            "and which have zero-derivative." << endl;
+        cout << "Try to decrease aDelta variable's value. ---" << endl;
+#endif
+        return Standard_False;
+        }
+        
+      F1 = FF1(2);
+      F2 = FF2(2);
+
+      Df(2,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU);
+      }//else of if(myU-myUinfium < 2*aDeltaU)
+/////////////////////////// Search of DF2_u derivative (end) ///////////////////
+    
+    //Return to previous values
+    myU = myU_old;
+    myV = myV_old;
+    myP1 = myP1_old;
+    myP2 = myP2_old;
+    myDu = myDu_old;
+    myDv = myDv_old;
+    }//if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
+  else
+    {
+    Standard_Real Ndv = myDv.Magnitude();
+    Df(2,2) = Ndv + (P1P2.Dot(Dvv)/Ndv) - F(2)*(myDv.Dot(Dvv)/(Ndv*Ndv));
+    Df(2,1) = -myDu.Dot(myDv)/Ndv;    
+    }//else of if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
   
-  Df(1,1) = - Ndu + (P1P2.Dot(Duu)/Ndu) - F(1)*(myDu.Dot(Duu)/(Ndu*Ndu));
-  Df(1,2) = myDv.Dot(myDu)/Ndu;
-  Df(2,1) = -myDu.Dot(myDv)/Ndv;
-  Df(2,2) = Ndv + (P1P2.Dot(Dvv)/Ndv) - F(2)*(myDv.Dot(Dvv)/(Ndv*Ndv));
   return Standard_True;
 
-}
+  }//end of function
 //=============================================================================
 
 Standard_Integer Extrema_FuncExtCC::GetStateNumber ()
@@ -166,11 +730,11 @@ Standard_Integer Extrema_FuncExtCC::GetStateNumber ()
   Vec P1P2 (myP1, myP2);
 
   Standard_Real mod = Du.Magnitude();
-  if(mod > Tol) {
+  if(mod > myTolC1) {
     Du /= mod;
   }
   mod = Dv.Magnitude();
-  if(mod > Tol) {
+  if(mod > myTolC2) {
     Dv /= mod;
   }
 
@@ -192,9 +756,12 @@ void Extrema_FuncExtCC::Points (const Standard_Integer N,
 }
 //=============================================================================
 
-
-
-
-
-
-
+void Extrema_FuncExtCC::SubIntervalInitialize(const math_Vector& theInfBound,
+                                              const math_Vector& theSupBound)
+  {
+  myUinfium = theInfBound(1);
+  myUsupremum = theSupBound(1);
+  myVinfium = theInfBound(2);
+  myVsupremum = theSupBound(2);
+  }
+//=============================================================================
\ No newline at end of file
index f66dacc..1ff093e 100755 (executable)
 // and conditions governing the rights and limitations under the License.
 
 //=============================================================================
-
-inline void Extrema_FuncExtCC::SetCurve (const Standard_Integer theRank, const Curve1& C)
-{
-  Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_FuncExtCC::SetCurve()")
-  if (theRank == 1) {myC1 = (Standard_Address)&C;}
-  else {myC2 = (Standard_Address)&C;}
-}
-
-//=============================================================================
-
 inline void Extrema_FuncExtCC::SetTolerance (const Standard_Real theTol)
 {
   myTol = theTol;
index dee78f8..6a7f9be 100755 (executable)
@@ -94,6 +94,14 @@ is
                 TypeMismatch from Standard;
                -- if N < 1 or N > NbExt(me).
 
+    SubIntervalInitialize(me: in out; theUfirst, theUlast: Real from Standard);
+       ---Purpose: Determines boundaries of subinterval for find of root.
+    
+    SearchOfTolerance(me: in out) returns Real from Standard;
+       ---Purpose: Computes a Tol value. If 1st derivative of curve
+       --          |D1|<Tol, it is considered D1=0.      
+    
+
 fields
     myP    : Pnt;
     myC    : Address from Standard;
@@ -108,5 +116,14 @@ fields
     myPinit: Boolean;
     myCinit: Boolean;
     myD1Init: Boolean;
+    
+    myTol: Real from Standard;  -- toolerance for derivate
+
+--Supremum of search 1st non-zero derivative    
+    myMaxDerivOrder: Integer from Standard;
+    
+--boundaries of subinterval for find of root
+    myUinfium, myUsupremum: Real from Standard;
+    
 
 end FuncExtPC;
index 7f4bc70..9869e23 100755 (executable)
 
 #include <Standard_TypeMismatch.hxx>
 #include <Precision.hxx>
-#define delta 1.e-9
-#define Tol   1.e-20
+static const Standard_Real TolFactor = 1.e-12;
+static const Standard_Real MinTol    = 1.e-20;
+static const Standard_Real MinStep   = 1e-7;
+
+static const Standard_Integer MaxOrder = 3;
+
+
 
 /*-----------------------------------------------------------------------------
  Fonction permettant de rechercher une distance extremale entre un point P et
@@ -36,6 +41,35 @@ les algorithmes math_FunctionRoot et math_FunctionRoots.
           = ||D1c(u)|| ** 2 + (C(u)-P).D2c(u) }
 ----------------------------------------------------------------------------*/
 
+Standard_Real Extrema_FuncExtPC::SearchOfTolerance()
+  {
+  const Standard_Integer NPoint = 10;
+  const Standard_Real aStep = (myUsupremum - myUinfium)/(Standard_Real)NPoint;
+  
+  Standard_Integer aNum = 0;
+  Standard_Real aMax = -Precision::Infinite();  //Maximum value of 1st derivative
+                                                //(it is computed with using NPoint point)
+    
+  do
+    {
+    Standard_Real u = myUinfium + aNum*aStep;    //parameter for every point
+    if(u > myUsupremum)
+      u = myUsupremum;
+    
+    Pnt Ptemp;  //empty point (is not used below)
+    Vec VDer;   // 1st derivative vector
+    Tool::D1(*((Curve*)myC), u, Ptemp, VDer);
+    Standard_Real vm = VDer.Magnitude();
+    if(vm > aMax)
+      aMax = vm;
+    }
+  while(++aNum < NPoint+1);
+  
+  return Max(aMax*TolFactor,MinTol);
+  
+  }
+
+//=============================================================================
 
 Extrema_FuncExtPC::Extrema_FuncExtPC():
 myU(0.),
@@ -44,31 +78,67 @@ myD1f(0.)
   myPinit = Standard_False;
   myCinit = Standard_False;
   myD1Init = Standard_False;
+  
+  SubIntervalInitialize(0.0,0.0);
+  myMaxDerivOrder = 0;
+  myTol=MinTol;
 }
 
 //=============================================================================
 
 Extrema_FuncExtPC::Extrema_FuncExtPC (const Pnt& P, 
-                                     const Curve& C):
-myU(0.),
-myD1f(0.) 
-{
+                                     const Curve& C): myU(0.), myD1f(0.)
+       {
   myP = P;
   myC = (Standard_Address)&C;
   myPinit = Standard_True;
   myCinit = Standard_True;
   myD1Init = Standard_False;
-}
+  
+  SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)),
+                             Tool::LastParameter(*((Curve*)myC)));
+                             
+  switch(Tool::GetType(*((Curve*)myC)))
+    {
+    case GeomAbs_BezierCurve:
+    case GeomAbs_BSplineCurve:
+    case GeomAbs_OtherCurve:
+      myMaxDerivOrder = MaxOrder;
+      myTol = SearchOfTolerance();
+      break;
+    default:
+      myMaxDerivOrder = 0;
+      myTol=MinTol;
+      break;
+    }
+  }
 //=============================================================================
 
 void Extrema_FuncExtPC::Initialize(const Curve& C)
-{
+  {
   myC = (Standard_Address)&C;
   myCinit = Standard_True;
   myPoint.Clear();
   mySqDist.Clear();
   myIsMin.Clear();
-}
+  
+  SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)),
+                             Tool::LastParameter(*((Curve*)myC)));
+                             
+  switch(Tool::GetType(*((Curve*)myC)))
+    {
+    case GeomAbs_BezierCurve:
+    case GeomAbs_BSplineCurve:
+    case GeomAbs_OtherCurve:
+      myMaxDerivOrder = MaxOrder;
+      myTol = SearchOfTolerance();
+      break;
+    default:
+      myMaxDerivOrder = 0;
+      myTol=MinTol;
+      break;
+    }
+  }
 
 //=============================================================================
 
@@ -83,30 +153,106 @@ void Extrema_FuncExtPC::SetPoint(const Pnt& P)
 
 //=============================================================================
 
+
 Standard_Boolean Extrema_FuncExtPC::Value (const Standard_Real U, Standard_Real& F)
 {
-  if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
+  if (!myPinit || !myCinit)
+    Standard_TypeMismatch::Raise("No init");
+  
   myU = U;
   Vec D1c;
   Tool::D1(*((Curve*)myC),myU,myPc,D1c);
   Standard_Real Ndu = D1c.Magnitude();
-  if (Ndu <= Tol) { // Cas Singulier (PMN 22/04/1998)
-       Pnt P1, P2;
-       P2 = Tool::Value(*((Curve*)myC),myU + delta);
-       P1 = Tool::Value(*((Curve*)myC),myU - delta);
-       Vec V(P1,P2);
-       D1c = V;
-       Ndu = D1c.Magnitude();
-       if (Ndu <= Tol) {
-         Vec aD2;
-         Tool::D2(*((Curve*)myC),myU,myPc,D1c,aD2);    
-         Ndu = aD2.Magnitude();
-        
-         if(Ndu  <= Tol)
-          return Standard_False;
-         D1c = aD2;    
+
+  if(myMaxDerivOrder != 0)
+    {
+    if (Ndu <= myTol) // Cas Singulier (PMN 22/04/1998)
+      {
+      const Standard_Real DivisionFactor = 1.e-3;
+      Standard_Real du;
+      if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
+        du = 0.0;
+      else
+        du = myUsupremum-myUinfium;
+        
+      const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
+//Derivative is approximated by Taylor-series
+      
+      Standard_Integer n = 1; //Derivative order
+      Vec V;
+      Standard_Boolean IsDeriveFound;
+      
+      do
+        {
+        V = Tool::DN(*((Curve*)myC),myU,++n);
+        Ndu = V.Magnitude();
+        IsDeriveFound = (Ndu > myTol);
+        }
+      while(!IsDeriveFound && n < myMaxDerivOrder);
+      
+      if(IsDeriveFound)
+        {
+        Standard_Real u;
+        
+        if(myU-myUinfium < aDelta)
+          u = myU+aDelta;
+        else
+          u = myU-aDelta;
+        
+        Pnt P1, P2;
+        Tool::D0(*((Curve*)myC),Min(myU, u),P1);
+        Tool::D0(*((Curve*)myC),Max(myU, u),P2);
+        
+        Vec V1(P1,P2);
+        Standard_Real aDirFactor = V.Dot(V1);
+        
+        if(aDirFactor < 0.0)
+          D1c = -V;
+        else
+          D1c = V;
+        }//if(IsDeriveFound)
+      else
+        {
+//Derivative is approximated by three points
+
+        Pnt Ptemp; //(0,0,0)-coordinate
+        Pnt P1, P2, P3;
+        Standard_Boolean IsParameterGrown;
+                  
+        if(myU-myUinfium < 2*aDelta)
+          {
+          Tool::D0(*((Curve*)myC),myU,P1);
+          Tool::D0(*((Curve*)myC),myU+aDelta,P2);
+          Tool::D0(*((Curve*)myC),myU+2*aDelta,P3);
+          IsParameterGrown = Standard_True;
+          }
+        else
+          {
+          Tool::D0(*((Curve*)myC),myU-2*aDelta,P1);
+          Tool::D0(*((Curve*)myC),myU-aDelta,P2);
+          Tool::D0(*((Curve*)myC),myU,P3);
+          IsParameterGrown = Standard_False;
+          }
+          
+        Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
+        
+        if(IsParameterGrown)
+          D1c=-3*V1+4*V2-V3;
+        else
+          D1c=V1-4*V2+3*V3;
+        }
+        Ndu = D1c.Magnitude();
+      }//(if (Ndu <= myTol)) condition
+    }//if(myMaxDerivOrder != 0)
+
+  if (Ndu <= MinTol)
+    {
+#ifdef DEB
+    cout << "+++Function Extrema_FuncExtPC::Value(...)." << endl;
+    cout << "Warning: 1st derivative is equal to zero!---"<<endl;
+#endif
+    return Standard_False;
     }
-  }
   
   Vec PPc (myP,myPc);
   F = PPc.Dot(D1c)/Ndu;
@@ -125,36 +271,104 @@ Standard_Boolean Extrema_FuncExtPC::Derivative (const Standard_Real U, Standard_
 //=============================================================================
 
 Standard_Boolean Extrema_FuncExtPC::Values (const Standard_Real U, Standard_Real& F, Standard_Real& D1f)
-{
-  if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
+  {
+  if (!myPinit || !myCinit)
+    Standard_TypeMismatch::Raise("No init");
+    
+  Pnt myPc_old = myPc, myP_old = myP;
+
+  if(Value(U,F) == Standard_False)
+    {
+#ifdef DEB
+    cout << "+++Function Extrema_FuncExtPC::Values(...)." << endl;
+    cout << "Warning: No function value found!---"<<endl;
+#endif
+
+    myD1Init = Standard_False;
+    return Standard_False;
+    }
+    
   myU = U;
+  myPc = myPc_old;
+  myP = myP_old;
+    
   Vec D1c,D2c;
   Tool::D2(*((Curve*)myC),myU,myPc,D1c,D2c);
 
   Standard_Real Ndu = D1c.Magnitude();
-  if (Ndu <= Tol) {// Cas Singulier (PMN 22/04/1998)
-    Pnt P1, P2;
-    Vec V1;
-    Tool::D1(*((Curve*)myC),myU+delta, P2, V1);
-    Tool::D1(*((Curve*)myC),myU-delta, P1, D2c);
-    Vec V(P1,P2);
-    D1c = V;
-    D2c -= V1;
-    Ndu = D1c.Magnitude();
-    if (Ndu <= Tol) {
-      myD1Init = Standard_False;
-      return Standard_False;
+  if (Ndu <= myTol) // Cas Singulier (PMN 22/04/1998)
+    {
+//Derivative is approximated by three points
+
+//Attention:  aDelta value must be greater than same value for "Value(...)"
+//            function to avoid of points' collisions.
+    const Standard_Real DivisionFactor = 0.01;
+    Standard_Real du;
+    if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst())) 
+      du = 0.0;
+    else
+      du = myUsupremum-myUinfium;
+      
+    const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
+    
+    Standard_Real F1, F2, F3;
+    
+    if(myU-myUinfium < 2*aDelta)
+      {
+      F1=F;
+      const Standard_Real U1 = myU, U2 = myU+aDelta, U3 = myU+2*aDelta;
+      
+      if(!((Value(U2,F2)) && (Value(U3,F3))))
+        {
+#ifdef DEB
+        cout << "+++ Function Extrema_FuncExtPC::Values(...)" << endl;
+        cout << "There are many points close to singularity points "
+            "and which have zero-derivative." << endl;
+        cout << "Try to decrease aDelta variable's value. ---" << endl;
+#endif
+        myD1Init = Standard_False;
+        return Standard_False;
+        }
+      
+//After calling of Value(...) function variable myU will be redeterminated.
+//So we must return it previous value.
+      D1f=(-3*F1+4*F2-F3)/(2.0*aDelta);
+      }
+    else
+      {
+      F3 = F;
+      const Standard_Real U1 = myU-2*aDelta, U2 = myU-aDelta, U3 = myU;
+      
+      if(!((Value(U2,F2)) && (Value(U1,F1))))
+        {
+#ifdef DEB
+        cout << "+++ Function Extrema_FuncExtPC::Values(...)" << endl;
+        cout << "There are many points close to singularity points "
+            "and which have zero-derivative." << endl;
+        cout << "Try to decrease aDelta variable's value. ---" << endl;
+#endif
+        myD1Init = Standard_False;
+        return Standard_False;
+        }
+//After calling of Value(...) function variable myU will be redeterminated.
+//So we must return it previous value.
+      D1f=(F1-4*F2+3*F3)/(2.0*aDelta);
+      }
+      myU = U;
+      myPc = myPc_old;
+      myP = myP_old;
+    }
+  else
+    {
+    Vec PPc (myP,myPc);
+    D1f = Ndu + (PPc.Dot(D2c)/Ndu) - F*(D1c.Dot(D2c))/(Ndu*Ndu);
     }
-  }
-  
-  Vec PPc (myP,myPc);
-  F = PPc.Dot(D1c)/Ndu;
-  D1f = Ndu + (PPc.Dot(D2c)/Ndu) - F*(D1c.Dot(D2c))/(Ndu*Ndu);
 
   myD1f = D1f;
+
   myD1Init = Standard_True;
   return Standard_True;
-}
+  }
 //=============================================================================
 
 Standard_Integer Extrema_FuncExtPC::GetStateNumber ()
@@ -200,3 +414,8 @@ POnC Extrema_FuncExtPC::Point (const Standard_Integer N) const
 }
 //=============================================================================
 
+void Extrema_FuncExtPC::SubIntervalInitialize(const Standard_Real theUfirst, const Standard_Real theUlast)
+  {
+  myUinfium = theUfirst;
+  myUsupremum = theUlast;
+  }
\ No newline at end of file
index d7e8b84..cf21029 100755 (executable)
 #include <Precision.hxx>
 
 
+//=======================================================================
+//function :  IsParallelDot
+//purpose  :  This function returns True if angle between <theV1> and 
+//<theV2> vectors or between <theV1> and <-theV2> vectors is less than 
+//AngTol.
+//=======================================================================
+static Standard_Boolean IsParallelDot(  gp_Vec2d theV1,
+                                        gp_Vec2d theV2,
+                                        Standard_Real AngTol)
+  {
+  return Abs(theV1.Dot(theV2)) >= 
+            theV1.Magnitude()*theV2.Magnitude()*cos(AngTol);
+  }
+
 Extrema_GExtCC2d::Extrema_GExtCC2d() {}
 
 
@@ -535,7 +549,8 @@ void Extrema_GExtCC2d::Results(const Extrema_ECC2d& AlgExt,
          gp_Vec2d v1, v2;
          Tool1::D1(C1,U,p, v1);
          Tool2::D1(*((Curve2*)myC),U2,p, v2);
-         if (v1.IsParallel(v2, Precision::Angular())) {
+         if (IsParallelDot(v1, v2, Precision::Angular()))
+           {
            mynbext++;
            Val = AlgExt.SquareDistance(i);
            P1.SetValues(U, P1.Value());
@@ -543,7 +558,7 @@ void Extrema_GExtCC2d::Results(const Extrema_ECC2d& AlgExt,
            mySqDist.Append(Val);
            mypoints.Append(P1);
            mypoints.Append(P2);
-         }
+           }
 //  modified by NIZHNY-EAP Thu Jan 27 16:41:00 2000 ___END___
        }
       }
index ddaf85e..54b87d4 100755 (executable)
@@ -131,8 +131,9 @@ void Extrema_GExtPC::Perform(const ThePoint& P)
         if (myuinf >= anInfToCheck) anInfToCheck = myuinf;
         if (myusup <= aSupToCheck) aSupToCheck = myusup;
         if((aSupToCheck - anInfToCheck) <= mytolu) continue;
-
-        if (i != 1) {
+        
+        if (i != 1)
+          {
           TheCurveTool::D1(*((TheCurve*)myC), myintuinf, PP, V1);
           s1 = (TheVector(P, PP))*V1;
           if (s1*s2 < 0.0) {
@@ -145,6 +146,7 @@ void Extrema_GExtPC::Perform(const ThePoint& P)
           TheCurveTool::D1(*((TheCurve*)myC), myintusup, PP, V1);
           s2 = (TheVector(P, PP))*V1;
         }
+
         IntervalPerform(P);
         IntExtIsDone = IntExtIsDone || mydone;
       }
index cff9b1e..b24371c 100755 (executable)
@@ -184,6 +184,8 @@ b- Calcul des minima:
   UVinf(2) = aCache2->TrimFirstParameter();
   UVsup(1) = aCache1->TrimLastParameter();
   UVsup(2) = aCache2->TrimLastParameter();
+  
+  myF.SubIntervalInitialize(UVinf,UVsup);
 
 //     - des 'bords' du tableau TbDist2
   for (NoV = 0; NoV <= aNbV+1; NoV++) {
index 2f354a3..2d5a3b9 100755 (executable)
@@ -171,6 +171,7 @@ Methode:
 -----------------------------------------------------------------------------*/
 {
   myF.SetPoint(P);
+  myF.SubIntervalInitialize(myumin,myusup);
   myDone = Standard_False;
 
   math_FunctionRoots S (myF, myumin, myusup, mynbsample, mytolu, mytolF, mytolF);
index cec8bb1..77e3277 100755 (executable)
@@ -60,10 +60,11 @@ typedef gp_XYZ  XYZ;
 
 
 
+//ordre de derivation maximum pour la recherche de la premiere 
+//derivee non nulle
+static const int maxDerivOrder = 3;
+static const Standard_Real MinStep   = 1e-7;
 
-static const int MaxDegree = 9;
-   //ordre de derivation maximum pour la recherche de la premiere 
-   //derivee non nulle
 
 
 
@@ -297,7 +298,7 @@ void Geom_OffsetCurve::D2 (const Standard_Real U, Pnt& P, Vec& V1, Vec& V2) cons
 //purpose  : 
 //=======================================================================
 
-void Geom_OffsetCurve::D3 (const Standard_Real U, Pnt& P, Vec& V1, Vec& V2, Vec& V3) 
+void Geom_OffsetCurve::D3 (const Standard_Real theU, Pnt& P, Vec& theV1, Vec& V2, Vec& V3) 
 const {
 
 
@@ -314,22 +315,71 @@ const {
    //         (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir +
    //         (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir
 
+  const Standard_Real aTol = gp::Resolution();
+
+  Standard_Boolean IsDirectionChange = Standard_False;
+
+  basisCurve->D3 (theU, P, theV1, V2, V3);
+  Vec V4 = basisCurve->DN (theU, 4);
+  if(theV1.Magnitude() <= aTol)
+    {
+    const Standard_Real anUinfium   = basisCurve->FirstParameter();
+    const Standard_Real anUsupremum = basisCurve->LastParameter();
+
+    const Standard_Real DivisionFactor = 1.e-3;
+    Standard_Real du;
+    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) 
+      du = 0.0;
+    else
+      du = anUsupremum-anUinfium;
+      
+    const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
+//Derivative is approximated by Taylor-series
+    
+    Standard_Integer anIndex = 1; //Derivative order
+    Vec V;
+    
+    do
+      {
+      V =  basisCurve->DN(theU,++anIndex);
+      }
+    while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder);
+    
+    Standard_Real u;
+    
+    if(theU-anUinfium < aDelta)
+      u = theU+aDelta;
+    else
+      u = theU-aDelta;
+    
+    Pnt P1, P2;
+    basisCurve->D0(Min(theU, u),P1);
+    basisCurve->D0(Max(theU, u),P2);
+    
+    Vec V1(P1,P2);
+    Standard_Real aDirFactor = V.Dot(V1);
+    
+    if(aDirFactor < 0.0)
+      {
+      theV1 = -V;
+      V2 = -basisCurve->DN (theU, anIndex + 1);
+      V3 = -basisCurve->DN (theU, anIndex + 2);
+      V4 = -basisCurve->DN (theU, anIndex + 3);
+
+      IsDirectionChange = Standard_True;
+      }
+    else
+      {
+      theV1 = V;
+      V2 = basisCurve->DN (theU, anIndex + 1);
+      V3 = basisCurve->DN (theU, anIndex + 2);
+      V4 = basisCurve->DN (theU, anIndex + 3);
+      }
+    }//if(V1.Magnitude() <= aTol)
 
 
-  basisCurve->D3 (U, P, V1, V2, V3);
-  Vec V4 = basisCurve->DN (U, 4);
-  Standard_Integer Index = 2;
-  while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
-    V1 = basisCurve->DN (U, Index);
-    Index++;
-  }
-  if (Index != 2) {
-    V2 = basisCurve->DN (U, Index);
-    V3 = basisCurve->DN (U, Index + 1);
-    V4 = basisCurve->DN (U, Index + 2);
-  }
   XYZ OffsetDir = direction.XYZ();
-  XYZ Ndir      = (V1.XYZ()).Crossed (OffsetDir);
+  XYZ Ndir      = (theV1.XYZ()).Crossed (OffsetDir);
   XYZ DNdir     = (V2.XYZ()).Crossed (OffsetDir);
   XYZ D2Ndir    = (V3.XYZ()).Crossed (OffsetDir);
   XYZ D3Ndir    = (V4.XYZ()).Crossed (OffsetDir);
@@ -351,7 +401,12 @@ const {
     D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 -
                                  15.0*Dr*Dr*Dr/R6 - D3r));
     D3Ndir.Multiply (offsetValue/R);
+
+    if(IsDirectionChange)
+      V3=-V3;
+
     V3.Add (Vec(D3Ndir));
+
     // V2 = P" (U) :
     Standard_Real R4 = R2 * R2;
     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
@@ -362,7 +417,7 @@ const {
     DNdir.Multiply(R);
     DNdir.Subtract (Ndir.Multiplied (Dr/R));
     DNdir.Multiply (offsetValue/R2);
-    V1.Add (Vec(DNdir));
+    theV1.Add (Vec(DNdir));
   }
   else {
     // V3 = P"' (U) :
@@ -372,7 +427,12 @@ const {
     D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 
                 15.0*Dr*Dr*Dr/R7 - D3r));
     D3Ndir.Multiply (offsetValue);
+    
+    if(IsDirectionChange)
+      V3=-V3;
+
     V3.Add (Vec(D3Ndir));
+    
     // V2 = P" (U) :
     D2Ndir.Divide (R);
     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R3));
@@ -382,12 +442,10 @@ const {
     // V1 = P' (U) :
     DNdir.Multiply (offsetValue/R);
     DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));        
-    V1.Add (Vec(DNdir));
+    theV1.Add (Vec(DNdir));
   }
   //P (U) :
-  Ndir.Multiply (offsetValue/R);
-  Ndir.Add (P.XYZ());
-  P.SetXYZ (Ndir);
+  D0(theU,P);
 }
 
 
@@ -396,134 +454,30 @@ const {
 //purpose  : 
 //=======================================================================
 
-Vec Geom_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const {
-
+Vec Geom_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const 
+  {
+  Standard_RangeError_Raise_if (N < 1, "Exception: "
+                              "Geom_OffsetCurve::DN(...). N<1.");
 
-  if (N < 1)                  Standard_RangeError::Raise();
-  XYZ OffsetDir = direction.XYZ();
-  Pnt P;
-  Vec V1, V2, dummy;
-  if (N == 1) {
-    basisCurve->D2 (U, P, V1, V2);
-    Standard_Integer Index = 2;
-    while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
-      V1 = basisCurve->DN (U, Index);
-      Index++;
-    }
-    if (Index != 2)  V2 = basisCurve->DN (U, Index);
-    XYZ Ndir  = (V1.XYZ()).Crossed (OffsetDir);
-    XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
-    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())  Geom_UndefinedDerivative::Raise();
-      Ndir.Multiply (Dr/R);
-      DNdir.Multiply(R);
-      DNdir.Subtract (Ndir);
-      DNdir.Multiply (offsetValue/R2);
-      V1.Add (Vec(DNdir));
-    }
-    else {
-      Ndir.Multiply (offsetValue * Dr / R3);
-      DNdir.Multiply (offsetValue/R);
-      DNdir.Subtract (Ndir);        
-      V1.Add (Vec(DNdir));
-    }
-    dummy = V1;
-  }
-  else if (N == 2) {
-    Vec V3;
-    basisCurve->D3 (U, P, V1, V2, V3);
-    Standard_Integer Index = 2;
-    while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
-      V1 = basisCurve->DN (U, Index);
-      Index++;
-    }
-    if (Index != 2) {
-      V2 = basisCurve->DN (U, Index);
-      V3 = basisCurve->DN (U, Index + 1);
-    }
-    XYZ Ndir   = (V1.XYZ()).Crossed (OffsetDir);
-    XYZ DNdir  = (V2.XYZ()).Crossed (OffsetDir);
-    XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir);
-    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())  Geom_UndefinedDerivative::Raise();
-      Ndir.Multiply ((3.0 * Dr * Dr / R4) - (D2r/R2));
-      DNdir.Multiply (2.0 * Dr / R2);
-      D2Ndir.Subtract (DNdir);
-      D2Ndir.Subtract (Ndir);
-      D2Ndir.Multiply (offsetValue / R);
-      V2.Add (Vec(D2Ndir));
-    }
-    else {
-      Ndir.Multiply ((3.0 * Dr * Dr / R4) - (D2r / R2));
-      DNdir.Multiply (2.0 * Dr / R2);
-      D2Ndir.Divide (R);
-      D2Ndir.Subtract (DNdir);
-      D2Ndir.Subtract (Ndir);
-      D2Ndir.Multiply (offsetValue);
-      V2.Add (Vec(D2Ndir));
-    }
-    dummy = V2;
-  }
-  else if (N == 3) {
-    Vec V3;
-    basisCurve->D3 (U, P, V1, V2, V3);
-    Vec V4 = basisCurve->DN (U, 4);
-    Standard_Integer Index = 2;
-    while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
-      V1 = basisCurve->DN (U, Index);
-      Index++;
-    }
-    if (Index != 2) {
-      V2 = basisCurve->DN (U, Index);
-      V3 = basisCurve->DN (U, Index + 1);
-      V4 = basisCurve->DN (U, Index + 2);
-    }
-    XYZ Ndir = (V1.XYZ()).Crossed (OffsetDir);
-    XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
-    XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir);
-    XYZ D3Ndir = (V4.XYZ()).Crossed (OffsetDir);
-    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()) Geom_UndefinedDerivative::Raise();
-      D2Ndir.Multiply (3.0 * Dr / R2);
-      DNdir.Multiply (3.0 * ((D2r/R2) + (Dr*Dr)/R4));
-      Ndir.Multiply (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r);
-      D3Ndir.Subtract (D2Ndir);
-      D3Ndir.Subtract (DNdir);
-      D3Ndir.Add (Ndir);
-      D3Ndir.Multiply (offsetValue/R);
-      V3.Add (Vec(D3Ndir));
-    }
-    else {
-      D2Ndir.Multiply (3.0 * Dr / R3);
-      DNdir.Multiplied (3.0 * ((D2r/R3) + (Dr*Dr/R5)));
-      Ndir.Multiply (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r);
-      D3Ndir.Divide (R);
-      D3Ndir.Subtract (D2Ndir);
-      D3Ndir.Subtract (DNdir);
-      D3Ndir.Add (Ndir);
-      D3Ndir.Multiply (offsetValue);
-      V3.Add (Vec(D3Ndir));
-    }
-    dummy = V3;
+  gp_Vec VN, Vtemp;
+  gp_Pnt Ptemp;
+  switch (N)
+    {
+    case 1:
+      D1( U, Ptemp, VN);
+      break;
+    case 2:
+      D2( U, Ptemp, Vtemp, VN);
+      break;
+    case 3:
+      D3( U, Ptemp, Vtemp, Vtemp, VN);
+      break;
+    default:
+      Standard_NotImplemented::Raise("Exception: "
+        "Derivative order is greater than 3. Cannot compute of derivative.");
   }
-  else { Standard_NotImplemented::Raise();  }
-
-  return dummy;
+  
+  return VN;
 }
 
 //=======================================================================
@@ -531,24 +485,70 @@ Vec Geom_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const
 //purpose  : 
 //=======================================================================
 
-void  Geom_OffsetCurve::D0(const Standard_Real U,
-                           gp_Pnt& P,
-                           gp_Pnt& Pbasis,
-                           gp_Vec& V1basis)const 
-{
+void  Geom_OffsetCurve::D0(const Standard_Real theU, gp_Pnt& theP,
+                           gp_Pnt& thePbasis, gp_Vec& theV1basis)const 
+  {
+  const Standard_Real aTol = gp::Resolution();
 
-  basisCurve->D1 (U, Pbasis, V1basis);
-  Standard_Integer Index = 2;
-  while (V1basis.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
-    V1basis = basisCurve->DN (U, Index);
-    Index++;
-  }
-  XYZ Ndir = (V1basis.XYZ()).Crossed (direction.XYZ());
+  basisCurve->D1 (theU, thePbasis, theV1basis);
+  Standard_Real Ndu = theV1basis.Magnitude();
+  
+  if(Ndu <= aTol)
+    {
+    const Standard_Real anUinfium   = basisCurve->FirstParameter();
+    const Standard_Real anUsupremum = basisCurve->LastParameter();
+
+    const Standard_Real DivisionFactor = 1.e-3;
+    Standard_Real du;
+    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) 
+      du = 0.0;
+    else
+      du = anUsupremum-anUinfium;
+      
+    const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
+//Derivative is approximated by Taylor-series
+    
+    Standard_Integer anIndex = 1; //Derivative order
+    gp_Vec V;
+    
+    do
+      {
+      V =  basisCurve->DN(theU,++anIndex);
+      Ndu = V.Magnitude();
+      }
+    while((Ndu <= aTol) && anIndex < maxDerivOrder);
+    
+    Standard_Real u;
+    
+    if(theU-anUinfium < aDelta)
+      u = theU+aDelta;
+    else
+      u = theU-aDelta;
+    
+    gp_Pnt P1, P2;
+    basisCurve->D0(Min(theU, u),P1);
+    basisCurve->D0(Max(theU, u),P2);
+    
+    gp_Vec V1(P1,P2);
+    Standard_Real aDirFactor = V.Dot(V1);
+    
+    if(aDirFactor < 0.0)
+      theV1basis = -V;
+    else
+      theV1basis = V;
+
+    Ndu = theV1basis.Magnitude();
+    }//if(Ndu <= aTol)
+  
+  XYZ Ndir = (theV1basis.XYZ()).Crossed (direction.XYZ());
   Standard_Real R = Ndir.Modulus();
-  if (R <= gp::Resolution())  Geom_UndefinedValue::Raise();
+  if (R <= gp::Resolution())
+    Geom_UndefinedValue::Raise("Exception: Undefined normal vector "
+          "because tangent vector has zero-magnitude!");
+
   Ndir.Multiply (offsetValue/R);
-  Ndir.Add (Pbasis.XYZ());
-  P.SetXYZ(Ndir);
+  Ndir.Add (thePbasis.XYZ());
+  theP.SetXYZ(Ndir);
 }
 
 //=======================================================================
@@ -556,26 +556,76 @@ void  Geom_OffsetCurve::D0(const Standard_Real U,
 //purpose  : 
 //=======================================================================
 
-void Geom_OffsetCurve::D1 ( const Standard_Real U, 
+void Geom_OffsetCurve::D1 ( const Standard_Real theU, 
                             Pnt& P , Pnt& PBasis ,
-                            Vec& V1, Vec& V1basis, Vec& V2basis) const {
+                            Vec& theV1, Vec& V1basis, Vec& V2basis) const {
 
    // P(u) = p(u) + Offset * Ndir / R
    // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction)
 
    // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R -  Ndir * (DR/R))
 
-  basisCurve->D2 (U, PBasis, V1basis, V2basis);
-  V1 = V1basis;
+  const Standard_Real aTol = gp::Resolution();
+
+  basisCurve->D2 (theU, PBasis, V1basis, V2basis);
+  theV1 = V1basis;
   Vec V2 = V2basis;
-  Standard_Integer Index = 2;
-  while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
-    V1 = basisCurve->DN (U, Index);
-    Index++;
-  }
-  if (Index != 2)  V2 = basisCurve->DN (U, Index);
+
+  if(theV1.Magnitude() <= aTol)
+    {
+    const Standard_Real anUinfium   = basisCurve->FirstParameter();
+    const Standard_Real anUsupremum = basisCurve->LastParameter();
+
+    const Standard_Real DivisionFactor = 1.e-3;
+    Standard_Real du;
+    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) 
+      du = 0.0;
+    else
+      du = anUsupremum-anUinfium;
+      
+    const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
+//Derivative is approximated by Taylor-series
+    
+    Standard_Integer anIndex = 1; //Derivative order
+    Vec V;
+    
+    do
+      {
+      V =  basisCurve->DN(theU,++anIndex);
+      }
+    while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder);
+    
+    Standard_Real u;
+    
+    if(theU-anUinfium < aDelta)
+      u = theU+aDelta;
+    else
+      u = theU-aDelta;
+    
+    Pnt P1, P2;
+    basisCurve->D0(Min(theU, u),P1);
+    basisCurve->D0(Max(theU, u),P2);
+    
+    Vec V1(P1,P2);
+    Standard_Real aDirFactor = V.Dot(V1);
+    
+    if(aDirFactor < 0.0)
+      {
+      theV1 = -V;
+      V2 = - basisCurve->DN (theU, anIndex+1);
+      }
+    else
+      {
+      theV1 = V;
+      V2 = basisCurve->DN (theU, anIndex+1);
+      }
+    
+    V2basis = V2;
+    V1basis = theV1;
+    }//if(theV1.Magnitude() <= aTol)
+
   XYZ OffsetDir = direction.XYZ();
-  XYZ Ndir  = (V1.XYZ()).Crossed (OffsetDir);
+  XYZ Ndir  = (theV1.XYZ()).Crossed (OffsetDir);
   XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir);
   Standard_Real R2 = Ndir.SquareModulus();
   Standard_Real R  = Sqrt (R2);
@@ -587,18 +637,16 @@ void Geom_OffsetCurve::D1 ( const Standard_Real U,
     DNdir.Multiply(R);
     DNdir.Subtract (Ndir.Multiplied (Dr/R));
     DNdir.Multiply (offsetValue/R2);
-    V1.Add (Vec(DNdir));
+    theV1.Add (Vec(DNdir));
   }
   else {
     // Same computation as IICURV in EUCLID-IS because the stability is
     // better
     DNdir.Multiply (offsetValue/R);
     DNdir.Subtract (Ndir.Multiplied (offsetValue * Dr/R3));        
-    V1.Add (Vec(DNdir));
+    theV1.Add (Vec(DNdir));
   }
-  Ndir.Multiply (offsetValue/R);
-  Ndir.Add (PBasis.XYZ());
-  P.SetXYZ (Ndir);
+  D0(theU,P);
 }
 
 
@@ -607,9 +655,9 @@ void Geom_OffsetCurve::D1 ( const Standard_Real U,
 //purpose  : 
 //=======================================================================
 
-void Geom_OffsetCurve::D2 (const Standard_Real U, 
+void Geom_OffsetCurve::D2 (const Standard_Real theU, 
                            Pnt& P      , Pnt& PBasis ,
-                           Vec& V1     , Vec& V2     , 
+                           Vec& theV1     , Vec& V2     , 
                            Vec& V1basis, Vec& V2basis, Vec& V3basis) const {
 
    // P(u) = p(u) + Offset * Ndir / R
@@ -620,21 +668,75 @@ void Geom_OffsetCurve::D2 (const Standard_Real U,
    // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) +
    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 
-  basisCurve->D3 (U, PBasis, V1basis, V2basis, V3basis);
-  Standard_Integer Index = 2;
-  V1     = V1basis;
+  const Standard_Real aTol = gp::Resolution();
+
+  Standard_Boolean IsDirectionChange = Standard_False;
+
+  basisCurve->D3 (theU, PBasis, V1basis, V2basis, V3basis);
+
+  theV1  = V1basis;
   V2     = V2basis;
   Vec V3 = V3basis;
-  while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
-    V1 = basisCurve->DN (U, Index);
-    Index++;
-  }
-  if (Index != 2) {
-    V2 = basisCurve->DN (U, Index);
-    V3 = basisCurve->DN (U, Index + 1);
-  }
+  
+  if(theV1.Magnitude() <= aTol)
+    {
+    const Standard_Real anUinfium   = basisCurve->FirstParameter();
+    const Standard_Real anUsupremum = basisCurve->LastParameter();
+
+    const Standard_Real DivisionFactor = 1.e-3;
+    Standard_Real du;
+    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) 
+      du = 0.0;
+    else
+      du = anUsupremum-anUinfium;
+      
+    const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
+//Derivative is approximated by Taylor-series
+    
+    Standard_Integer anIndex = 1; //Derivative order
+    Vec V;
+    
+    do
+      {
+      V =  basisCurve->DN(theU,++anIndex);
+      }
+    while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder);
+    
+    Standard_Real u;
+    
+    if(theU-anUinfium < aDelta)
+      u = theU+aDelta;
+    else
+      u = theU-aDelta;
+    
+    Pnt P1, P2;
+    basisCurve->D0(Min(theU, u),P1);
+    basisCurve->D0(Max(theU, u),P2);
+    
+    Vec V1(P1,P2);
+    Standard_Real aDirFactor = V.Dot(V1);
+    
+    if(aDirFactor < 0.0)
+      {
+      theV1 = -V;
+      V2 = -basisCurve->DN (theU, anIndex+1);
+      V3 = -basisCurve->DN (theU, anIndex + 2);
+
+      IsDirectionChange = Standard_True;
+      }
+    else
+      {
+      theV1 = V;
+      V2 = basisCurve->DN (theU, anIndex+1);
+      V3 = basisCurve->DN (theU, anIndex + 2);
+      }
+    
+    V2basis = V2;
+    V1basis = theV1;
+    }//if(V1.Magnitude() <= aTol)
+  
   XYZ OffsetDir = direction.XYZ();
-  XYZ Ndir   = (V1.XYZ()).Crossed (OffsetDir);
+  XYZ Ndir   = (theV1.XYZ()).Crossed (OffsetDir);
   XYZ DNdir  = (V2.XYZ()).Crossed (OffsetDir);
   XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir);
   Standard_Real R2    = Ndir.SquareModulus();
@@ -644,6 +746,7 @@ void Geom_OffsetCurve::D2 (const Standard_Real U,
   Standard_Real R5    = R3 * R2;
   Standard_Real Dr    = Ndir.Dot (DNdir);
   Standard_Real D2r   = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir);
+  
   if (R5 <= gp::Resolution()) {
     //We try another computation but the stability is not very good
     //dixit ISG.
@@ -653,12 +756,17 @@ void Geom_OffsetCurve::D2 (const Standard_Real U,
     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
     D2Ndir.Add (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
     D2Ndir.Multiply (offsetValue / R);
+    
+    if(IsDirectionChange)
+      V2=-V2;
+
     V2.Add (Vec(D2Ndir));
+        
     // V1 = P' (U) :
     DNdir.Multiply(R);
     DNdir.Subtract (Ndir.Multiplied (Dr/R));
     DNdir.Multiply (offsetValue/R2);
-    V1.Add (Vec(DNdir));
+    theV1.Add (Vec(DNdir));
   }
   else {
     // Same computation as IICURV in EUCLID-IS because the stability is
@@ -670,16 +778,19 @@ void Geom_OffsetCurve::D2 (const Standard_Real U,
                      offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))
                                      )
                     );
+    
+    if(IsDirectionChange)
+      V2=-V2;
+
     V2.Add (Vec(D2Ndir));
+    
     // V1 = P' (U) :
     DNdir.Multiply (offsetValue/R);
     DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));        
-    V1.Add (Vec(DNdir));
+    theV1.Add (Vec(DNdir));
   }
   //P (U) :
-  Ndir.Multiply (offsetValue/R);
-  Ndir.Add (PBasis.XYZ());
-  P.SetXYZ (Ndir);
+  D0(theU,P);
 }
 
 
@@ -717,23 +828,15 @@ Standard_Real Geom_OffsetCurve::Offset () const { return offsetValue; }
 //purpose  : 
 //=======================================================================
 
-void Geom_OffsetCurve::Value (
-const Standard_Real U, Pnt& P, Pnt& PBasis,  Vec& V1basis) const {
+void Geom_OffsetCurve::Value (const Standard_Real theU, Pnt& theP, 
+                              Pnt& thePbasis,  Vec& theV1basis) const 
+  {
+  if (basisCurve->Continuity() == GeomAbs_C0)
+    Geom_UndefinedValue::Raise("Exception: Basis curve is C0 continuity!");
 
-  if (basisCurve->Continuity() == GeomAbs_C0)  Geom_UndefinedValue::Raise();
-  basisCurve->D1 (U, PBasis, V1basis);
-  Standard_Integer Index = 2;
-  while (V1basis.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
-    V1basis = basisCurve->DN (U, Index);
-    Index++;
+  basisCurve->D1(theU, thePbasis, theV1basis);
+  D0(theU,theP);
   }
-  XYZ Ndir = (V1basis.XYZ()).Crossed (direction.XYZ());
-  Standard_Real R = Ndir.Modulus();
-  if (R <= gp::Resolution())  Geom_UndefinedValue::Raise();
-  Ndir.Multiply (offsetValue/R);
-  Ndir.Add (PBasis.XYZ());
-  P.SetXYZ (Ndir);
-}
 
 
 //=======================================================================
index 583e88c..dfd5fc0 100755 (executable)
@@ -53,11 +53,10 @@ typedef gp_Trsf2d Trsf2d;
 typedef gp_XY     XY;
 
 
-
-static const int MaxDegree = 9;
 //ordre de derivation maximum pour la recherche de la premiere 
 //derivee non nulle
-
+static const int maxDerivOrder = 3;
+static const Standard_Real MinStep   = 1e-7;
 
 
 //=======================================================================
@@ -192,47 +191,141 @@ GeomAbs_Shape Geom2d_OffsetCurve::Continuity () const
 //purpose  : 
 //=======================================================================
 
-void Geom2d_OffsetCurve::D0 (const Standard_Real   U,
-                                  Pnt2d& P ) const 
-{
-  Vec2d V1;
+void Geom2d_OffsetCurve::D0 (const Standard_Real   theU,
+                                  Pnt2d& theP ) const 
+  {
+  const Standard_Real aTol = gp::Resolution();
 
-  basisCurve->D1 (U, P, V1);
-  Standard_Integer Index = 2;
-  while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
-    V1 = basisCurve->DN (U, Index);
-    Index++;
+  Vec2d vD1;
+
+  basisCurve->D1 (theU, theP, vD1);
+  Standard_Real Ndu = vD1.Magnitude();
+
+  if(Ndu <= aTol)
+    {
+    const Standard_Real anUinfium   = basisCurve->FirstParameter();
+    const Standard_Real anUsupremum = basisCurve->LastParameter();
+
+    const Standard_Real DivisionFactor = 1.e-3;
+    Standard_Real du;
+    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) 
+      du = 0.0;
+    else
+      du = anUsupremum-anUinfium;
+      
+    const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
+//Derivative is approximated by Taylor-series
+    
+    Standard_Integer anIndex = 1; //Derivative order
+    Vec2d V;
+    
+    do
+      {
+      V =  basisCurve->DN(theU,++anIndex);
+      Ndu = V.Magnitude();
+      }
+    while((Ndu <= aTol) && anIndex < maxDerivOrder);
+    
+    Standard_Real u;
+    
+    if(theU-anUinfium < aDelta)
+      u = theU+aDelta;
+    else
+      u = theU-aDelta;
+    
+    Pnt2d P1, P2;
+    basisCurve->D0(Min(theU, u),P1);
+    basisCurve->D0(Max(theU, u),P2);
+    
+    Vec2d V1(P1,P2);
+    Standard_Real aDirFactor = V.Dot(V1);
+    
+    if(aDirFactor < 0.0)
+      vD1 = -V;
+    else
+      vD1 = V;
+
+    Ndu = vD1.Magnitude();
+    }//if(Ndu <= aTol)
+
+  if (Ndu <= aTol)
+    Geom2d_UndefinedValue::Raise("Exception: Undefined normal vector "
+          "because tangent vector has zero-magnitude!");
+  
+  Standard_Real A = vD1.Y();
+  Standard_Real B = - vD1.X();
+  A = A * offsetValue/Ndu;
+  B = B * offsetValue/Ndu;
+  theP.SetCoord(theP.X() + A, theP.Y() + B);
   }
-  Standard_Real A = V1.Y();
-  Standard_Real B = - V1.X();
-  Standard_Real R = Sqrt(A*A + B * B);
-  if (R <= gp::Resolution())  Geom2d_UndefinedValue::Raise();
-  A = A * offsetValue/R;
-  B = B * offsetValue/R;
-  P.SetCoord(P.X() + A, P.Y() + B);
-}
 
 //=======================================================================
 //function : D1
 //purpose  : 
 //=======================================================================
-
-void Geom2d_OffsetCurve::D1 (const Standard_Real U, Pnt2d& P, Vec2d& V1) const {
-
+void Geom2d_OffsetCurve::D1 (const Standard_Real theU, Pnt2d& P, Vec2d& theV1) 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))
 
+  const Standard_Real aTol = gp::Resolution();
+
   Vec2d V2;
-  basisCurve->D2 (U, P, V1, V2);
-  Standard_Integer Index = 2;
-  while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
-    V1 = basisCurve->DN (U, Index);
-    Index++;
-  }
-  if (Index != 2) { V2 = basisCurve->DN (U, Index); }
-  XY Ndir  (V1.Y(), -V1.X());
+  basisCurve->D2 (theU, P, theV1, V2);
+  
+  if(theV1.Magnitude() <= aTol)
+    {
+    const Standard_Real anUinfium   = basisCurve->FirstParameter();
+    const Standard_Real anUsupremum = basisCurve->LastParameter();
+
+    const Standard_Real DivisionFactor = 1.e-3;
+    Standard_Real du;
+    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) 
+      du = 0.0;
+    else
+      du = anUsupremum-anUinfium;
+      
+    const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
+//Derivative is approximated by Taylor-series
+    
+    Standard_Integer anIndex = 1; //Derivative order
+    Vec2d V;
+    
+    do
+      {
+      V =  basisCurve->DN(theU,++anIndex);
+      }
+    while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder);
+    
+    Standard_Real u;
+    
+    if(theU-anUinfium < aDelta)
+      u = theU+aDelta;
+    else
+      u = theU-aDelta;
+    
+    Pnt2d P1, P2;
+    basisCurve->D0(Min(theU, u),P1);
+    basisCurve->D0(Max(theU, u),P2);
+    
+    Vec2d V1(P1,P2);
+    Standard_Real aDirFactor = V.Dot(V1);
+    
+    if(aDirFactor < 0.0)
+      {
+      theV1 = -V;
+      V2 = - basisCurve->DN (theU, anIndex+1);
+      }
+    else
+      {
+      theV1 = V;
+      V2 = basisCurve->DN (theU, anIndex+1);
+      }
+    }//if(theV1.Magnitude() <= aTol)
+
+  XY Ndir  (theV1.Y(), -theV1.X());
   XY DNdir (V2.Y(), -V2.X());
   Standard_Real R2 = Ndir.SquareModulus();
   Standard_Real R  = Sqrt (R2);
@@ -244,18 +337,17 @@ void Geom2d_OffsetCurve::D1 (const Standard_Real U, Pnt2d& P, Vec2d& V1) const {
     DNdir.Multiply(R);
     DNdir.Subtract (Ndir.Multiplied (Dr/R));
     DNdir.Multiply (offsetValue/R2);
-    V1.Add (Vec2d(DNdir));
+    theV1.Add (Vec2d(DNdir));
   }
   else {
     // Same computation as IICURV in EUCLID-IS because the stability is
     // better
     DNdir.Multiply (offsetValue/R);
     DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));        
-    V1.Add (Vec2d(DNdir));
+    theV1.Add (Vec2d(DNdir));
   }
-  Ndir.Multiply (offsetValue/R);
-  Ndir.Add (P.XY());
-  P.SetXY (Ndir);
+
+  D0(theU, P);
 }
 
 //=======================================================================
@@ -263,9 +355,9 @@ void Geom2d_OffsetCurve::D1 (const Standard_Real U, Pnt2d& P, Vec2d& V1) const {
 //purpose  : 
 //=======================================================================
 
-void Geom2d_OffsetCurve::D2 (const Standard_Real U, 
+void Geom2d_OffsetCurve::D2 (const Standard_Real theU, 
                                   Pnt2d& P, 
-                                  Vec2d& V1, Vec2d& V2) const 
+                                  Vec2d& theV1, Vec2d& V2) const 
 {
    // P(u) = p(u) + Offset * Ndir / R
    // with R = || p' ^ Z|| and Ndir = P' ^ Z
@@ -276,17 +368,67 @@ void Geom2d_OffsetCurve::D2 (const Standard_Real U,
    //         Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2)))
 
   Vec2d V3;
-  basisCurve->D3 (U, P, V1, V2, V3);
-  Standard_Integer Index = 2;
-  while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
-    V1 = basisCurve->DN (U, Index);
-    Index++;
-  }
-  if (Index != 2) {
-    V2 = basisCurve->DN (U, Index);
-    V3 = basisCurve->DN (U, Index + 1);
-  }
-  XY Ndir (V1.Y(), -V1.X());
+  basisCurve->D3 (theU, P, theV1, V2, V3);
+
+  const Standard_Real aTol = gp::Resolution();
+
+  Standard_Boolean IsDirectionChange = Standard_False;
+
+  if(theV1.Magnitude() <= aTol)
+    {
+    const Standard_Real anUinfium   = basisCurve->FirstParameter();
+    const Standard_Real anUsupremum = basisCurve->LastParameter();
+
+    const Standard_Real DivisionFactor = 1.e-3;
+    Standard_Real du;
+    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) 
+      du = 0.0;
+    else
+      du = anUsupremum-anUinfium;
+      
+    const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
+//Derivative is approximated by Taylor-series
+    
+    Standard_Integer anIndex = 1; //Derivative order
+    Vec2d V;
+    
+    do
+      {
+      V =  basisCurve->DN(theU,++anIndex);
+      }
+    while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder);
+    
+    Standard_Real u;
+    
+    if(theU-anUinfium < aDelta)
+      u = theU+aDelta;
+    else
+      u = theU-aDelta;
+    
+    Pnt2d P1, P2;
+    basisCurve->D0(Min(theU, u),P1);
+    basisCurve->D0(Max(theU, u),P2);
+    
+    Vec2d V1(P1,P2);
+    Standard_Real aDirFactor = V.Dot(V1);
+    
+    if(aDirFactor < 0.0)
+      {
+      theV1 = -V;
+      V2 = -basisCurve->DN (theU, anIndex+1);
+      V3 = -basisCurve->DN (theU, anIndex + 2);
+
+      IsDirectionChange = Standard_True;
+      }
+    else
+      {
+      theV1 = V;
+      V2 = basisCurve->DN (theU, anIndex+1);
+      V3 = basisCurve->DN (theU, anIndex + 2);
+      }
+    }//if(V1.Magnitude() <= aTol)
+
+  XY Ndir (theV1.Y(), -theV1.X());
   XY DNdir (V2.Y(), -V2.X());
   XY D2Ndir (V3.Y(), -V3.X());
   Standard_Real R2  = Ndir.SquareModulus();
@@ -296,40 +438,54 @@ void Geom2d_OffsetCurve::D2 (const Standard_Real U,
   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 (R5 <= gp::Resolution())
+    {
     //We try another computation but the stability is not very good
     //dixit ISG.
-    if (R4 <= gp::Resolution()) { Geom2d_UndefinedDerivative::Raise(); }
+    if (R4 <= gp::Resolution())
+      {
+      Geom2d_UndefinedDerivative::Raise();
+      }
+    // V2 = P" (U) :
+    Standard_Real R4 = R2 * R2;
+    D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
+    D2Ndir.Add (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
+    D2Ndir.Multiply (offsetValue / R);
+
+    if(IsDirectionChange)
+      V2=-V2;
+
+    V2.Add (Vec2d(D2Ndir));
+
+    // V1 = P' (U) :
+    DNdir.Multiply(R);
+    DNdir.Subtract (Ndir.Multiplied (Dr/R));
+    DNdir.Multiply (offsetValue/R2);
+    theV1.Add (Vec2d(DNdir));
+    }
+  else
+    {
+    // Same computation as IICURV in EUCLID-IS because the stability is
+    // better.
     // V2 = P" (U) :
-     Standard_Real R4 = R2 * R2;
-     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
-     D2Ndir.Add (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
-     D2Ndir.Multiply (offsetValue / R);
-     V2.Add (Vec2d(D2Ndir));
-     // V1 = P' (U) :
-     DNdir.Multiply(R);
-     DNdir.Subtract (Ndir.Multiplied (Dr/R));
-     DNdir.Multiply (offsetValue/R2);
-     V1.Add (Vec2d(DNdir));
-   }
-   else {
-     // Same computation as IICURV in EUCLID-IS because the stability is
-     // better.
-     // V2 = P" (U) :
     D2Ndir.Multiply (offsetValue/R);
     D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
     D2Ndir.Add (Ndir.Multiplied 
                     (offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))));
+
+    if(IsDirectionChange)
+      V2=-V2;
+
     V2.Add (Vec2d(D2Ndir));
+
     // V1 = P' (U) 
-     DNdir.Multiply (offsetValue/R);
-     DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));        
-     V1.Add (Vec2d(DNdir));
-   }
-   //P (U) :
-   Ndir.Multiply (offsetValue/R);
-   Ndir.Add (P.XY());
-   P.SetXY (Ndir);
+    DNdir.Multiply (offsetValue/R);
+    DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));        
+    theV1.Add (Vec2d(DNdir));
+    }
+
+  //P (U) :
+  D0(theU, P);
 }
 
 
@@ -338,9 +494,9 @@ void Geom2d_OffsetCurve::D2 (const Standard_Real U,
 //purpose  : 
 //=======================================================================
 
-void Geom2d_OffsetCurve::D3 (const Standard_Real U, 
+void Geom2d_OffsetCurve::D3 (const Standard_Real theU, 
                                    Pnt2d& P, 
-                                   Vec2d& V1, Vec2d& V2, Vec2d& V3) const {
+                                   Vec2d& theV1, Vec2d& V2, Vec2d& V3) const {
 
 
    // P(u) = p(u) + Offset * Ndir / R
@@ -356,89 +512,149 @@ void Geom2d_OffsetCurve::D3 (const Standard_Real U,
    //         (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir +
    //         (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir
 
+  const Standard_Real aTol = gp::Resolution();
 
+  Standard_Boolean IsDirectionChange = Standard_False;
 
-   basisCurve->D3 (U, P, V1, V2, V3);
-   Vec2d V4 = basisCurve->DN (U, 4);
-   Standard_Integer Index = 2;
-   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
-     V1 = basisCurve->DN (U, Index);
-     Index++;
-   }
-   if (Index != 2) {
-     V2 = basisCurve->DN (U, Index);
-     V3 = basisCurve->DN (U, Index + 1);
-     V4 = basisCurve->DN (U, Index + 2);
-   }
-   XY Ndir   (V1.Y(), -V1.X());
-   XY DNdir  (V2.Y(), -V2.X());
-   XY D2Ndir (V3.Y(), -V3.X());
-   XY D3Ndir (V4.Y(), -V4.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()) {
-     //We try another computation but the stability is not very good
-     //dixit ISG.
-     if (R6 <= gp::Resolution()) Geom2d_UndefinedDerivative::Raise();
-     // V3 = P"' (U) :
-     D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * offsetValue * Dr / R2));
-     D3Ndir.Subtract (
-     (DNdir.Multiplied ((3.0 * offsetValue) * ((D2r/R2) + (Dr*Dr)/R4))));
-     D3Ndir.Add (Ndir.Multiplied (
-     (offsetValue * (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r))
-     ));
-     D3Ndir.Multiply (offsetValue/R);
-     V3.Add (Vec2d(D3Ndir));
-     // V2 = P" (U) :
-     Standard_Real R4 = R2 * R2;
-     D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
-     D2Ndir.Subtract (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
-     D2Ndir.Multiply (offsetValue / R);
-     V2.Add (Vec2d(D2Ndir));
-     // V1 = P' (U) :
-     DNdir.Multiply(R);
-     DNdir.Subtract (Ndir.Multiplied (Dr/R));
-     DNdir.Multiply (offsetValue/R2);
-     V1.Add (Vec2d(DNdir));
-   }
-   else {
-     // Same computation as IICURV in EUCLID-IS because the stability is
-     // better.
-     // V3 = P"' (U) :
-      D3Ndir.Multiply (offsetValue/R);
-     D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * offsetValue * Dr / R3));
-     D3Ndir.Subtract (DNdir.Multiplied (
-     ((3.0 * offsetValue) * ((D2r/R3) + (Dr*Dr)/R5))) );
-     D3Ndir.Add (Ndir.Multiplied (
-     (offsetValue * (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r))
-     ));
-     V3.Add (Vec2d(D3Ndir));
-     // V2 = P" (U) :
-     D2Ndir.Multiply (offsetValue/R);
-     D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
-     D2Ndir.Subtract (Ndir.Multiplied (
-                      offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))
-                                      )
-                     );
-     V2.Add (Vec2d(D2Ndir));
-     // V1 = P' (U) :
-     DNdir.Multiply (offsetValue/R);
-     DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));        
-     V1.Add (Vec2d(DNdir));
-   }
-   //P (U) :
-   Ndir.Multiply (offsetValue/R);
-   Ndir.Add (P.XY());
-   P.SetXY (Ndir);
-}
+  basisCurve->D3 (theU, P, theV1, V2, V3);
+  Vec2d V4 = basisCurve->DN (theU, 4);
+
+  if(theV1.Magnitude() <= aTol)
+    {
+    const Standard_Real anUinfium   = basisCurve->FirstParameter();
+    const Standard_Real anUsupremum = basisCurve->LastParameter();
+
+    const Standard_Real DivisionFactor = 1.e-3;
+    Standard_Real du;
+    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) 
+      du = 0.0;
+    else
+      du = anUsupremum-anUinfium;
+      
+    const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
+//Derivative is approximated by Taylor-series
+    
+    Standard_Integer anIndex = 1; //Derivative order
+    Vec2d V;
+    
+    do
+      {
+      V =  basisCurve->DN(theU,++anIndex);
+      }
+    while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder);
+    
+    Standard_Real u;
+    
+    if(theU-anUinfium < aDelta)
+      u = theU+aDelta;
+    else
+      u = theU-aDelta;
+    
+    Pnt2d P1, P2;
+    basisCurve->D0(Min(theU, u),P1);
+    basisCurve->D0(Max(theU, u),P2);
+    
+    Vec2d V1(P1,P2);
+    Standard_Real aDirFactor = V.Dot(V1);
+    
+    if(aDirFactor < 0.0)
+      {
+      theV1 = -V;
+      V2 = -basisCurve->DN (theU, anIndex + 1);
+      V3 = -basisCurve->DN (theU, anIndex + 2);
+      V4 = -basisCurve->DN (theU, anIndex + 3);
+
+      IsDirectionChange = Standard_True;
+      }
+    else
+      {
+      theV1 = V;
+      V2 = basisCurve->DN (theU, anIndex + 1);
+      V3 = basisCurve->DN (theU, anIndex + 2);
+      V4 = basisCurve->DN (theU, anIndex + 3);
+      }
+    }//if(V1.Magnitude() <= aTol)
+
+  XY Ndir   (theV1.Y(), -theV1.X());
+  XY DNdir  (V2.Y(), -V2.X());
+  XY D2Ndir (V3.Y(), -V3.X());
+  XY D3Ndir (V4.Y(), -V4.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()) 
+    {
+    //We try another computation but the stability is not very good
+    //dixit ISG.
+
+    if (R6 <= gp::Resolution())
+      Geom2d_UndefinedDerivative::Raise();
+
+    // V3 = P"' (U) :
+    D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * offsetValue * Dr / R2));
+    D3Ndir.Subtract (
+      (DNdir.Multiplied ((3.0 * offsetValue) * ((D2r/R2) + (Dr*Dr)/R4))));
+    D3Ndir.Add (Ndir.Multiplied (
+      (offsetValue * (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r))));
+    D3Ndir.Multiply (offsetValue/R);
+
+    if(IsDirectionChange)
+      V3=-V3;
+
+    V3.Add (Vec2d(D3Ndir));
+
+
+    // V2 = P" (U) :
+    Standard_Real R4 = R2 * R2;
+    D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2));
+    D2Ndir.Subtract (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2)));
+    D2Ndir.Multiply (offsetValue / R);
+    V2.Add (Vec2d(D2Ndir));
+    // V1 = P' (U) :
+    DNdir.Multiply(R);
+    DNdir.Subtract (Ndir.Multiplied (Dr/R));
+    DNdir.Multiply (offsetValue/R2);
+    theV1.Add (Vec2d(DNdir));
+    }
+  else
+    {
+    // Same computation as IICURV in EUCLID-IS because the stability is
+    // better.
+    // V3 = P"' (U) :
+    D3Ndir.Multiply (offsetValue/R);
+    D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * offsetValue * Dr / R3));
+    D3Ndir.Subtract (DNdir.Multiplied (
+      ((3.0 * offsetValue) * ((D2r/R3) + (Dr*Dr)/R5))) );
+    D3Ndir.Add (Ndir.Multiplied (
+      (offsetValue * (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r))));
+
+    if(IsDirectionChange)
+      V3=-V3;
+
+    V3.Add (Vec2d(D3Ndir));
+
+    // V2 = P" (U) :
+    D2Ndir.Multiply (offsetValue/R);
+    D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3));
+    D2Ndir.Subtract (Ndir.Multiplied (
+      offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3))));
+    V2.Add (Vec2d(D2Ndir));
+    // V1 = P' (U) :
+    DNdir.Multiply (offsetValue/R);
+    DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3));        
+    theV1.Add (Vec2d(DNdir));
+    }
+  //P (U) :
+  D0(theU, P);
+  }
 
 //=======================================================================
 //function : DN
@@ -448,7 +664,7 @@ void Geom2d_OffsetCurve::D3 (const Standard_Real U,
 Vec2d Geom2d_OffsetCurve::DN (const Standard_Real U, 
                               const Standard_Integer N) const 
 {
-  Standard_RangeError_Raise_if (N < 1, "Geom2d_OffsetCurve::DN()");
+Standard_RangeError_Raise_if (N < 1, "Exception: Geom2d_OffsetCurve::DN(). N<1.");
 
   gp_Vec2d VN, VBidon;
   gp_Pnt2d PBidon;
@@ -457,7 +673,8 @@ Vec2d Geom2d_OffsetCurve::DN (const Standard_Real U,
   case 2: D2( U, PBidon, VBidon, VN); break;
   case 3: D3( U, PBidon, VBidon, VBidon, VN); break;
   default:
-    Standard_NotImplemented::Raise();
+    Standard_NotImplemented::Raise("Exception: Derivative order is greater than 3. "
+      "Cannot compute of derivative.");
   }
   
   return VN;
@@ -469,25 +686,13 @@ Vec2d Geom2d_OffsetCurve::DN (const Standard_Real U,
 //purpose  : 
 //=======================================================================
 
-void Geom2d_OffsetCurve::Value (const Standard_Real U, 
-                               Pnt2d& P, Pnt2d& Pbasis,
-                               Vec2d& V1basis ) const 
-{
-
-  basisCurve->D1 (U, Pbasis, V1basis);
-  Standard_Integer Index = 2;
-  while (V1basis.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
-    V1basis = basisCurve->DN (U, Index);
-    Index++;
+void Geom2d_OffsetCurve::Value (const Standard_Real theU, 
+                                Pnt2d& theP, Pnt2d& thePbasis,
+                                Vec2d& theV1basis ) const 
+  {
+  basisCurve->D1(theU, thePbasis, theV1basis);
+  D0(theU,theP);
   }
-  Standard_Real A = V1basis.Y();
-  Standard_Real B = - V1basis.X();
-  Standard_Real R = Sqrt(A*A + B * B);
-  if (R <= gp::Resolution())  Geom2d_UndefinedValue::Raise();
-  A = A * offsetValue/R;
-  B = B * offsetValue/R;
-  P.SetCoord (A + Pbasis.X(), B + Pbasis.Y());
-}
 
 
 //=======================================================================
@@ -509,7 +714,7 @@ void Geom2d_OffsetCurve::D1 (const Standard_Real U,
    V1 = V1basis;
    Vec2d V2 = V2basis;
    Standard_Integer Index = 2;
-   while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
+   while (V1.Magnitude() <= gp::Resolution() && Index <= maxDerivOrder) {
      V1 = basisCurve->DN (U, Index);
      Index++;
    }
@@ -567,7 +772,7 @@ void Geom2d_OffsetCurve::D2 (const Standard_Real U,
   V1 = V1basis;
   V2 = V2basis;
   Vec2d V3 = V3basis;
-  while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) {
+  while (V1.Magnitude() <= gp::Resolution() && Index <= maxDerivOrder) {
     V1 = basisCurve->DN (U, Index);
     Index++;
   }
index 9c78ab2..5b4e8d4 100755 (executable)
@@ -44,7 +44,7 @@ uses
 is
 
 
-    TheType(myclass; C: IntCurveCurve)
+    GetType(myclass; C: IntCurveCurve)
        ---C++: inline
        returns CurveType from GeomAbs;
 
@@ -145,6 +145,15 @@ is
     D2 (myclass; C: IntCurveCurve; U: Real ; 
                  P: out Pnt2d; T,N: out Vec2d);
          ---C++: inline
+    
+    D3 (myclass; C: IntCurveCurve; U: Real ; 
+                 P: out Pnt2d; T,N,V: out Vec2d);
+         ---C++: inline
+    
+    DN(myclass; C: IntCurveCurve; U: Real ; 
+                 N: Integer from Standard)
+         returns Vec2d;
+         ---C++: inline
 
    NbIntervals(myclass ; C:  IntCurveCurve) 
         ---Purpose : output the number of interval of continuity C2 of 
@@ -165,6 +174,10 @@ is
         --           used if Type == Composite.
         ---C++: inline
 
+   Degree(myclass; C : IntCurveCurve) returns Integer from Standard;
+       ---C++: inline
+
+
 end CurveTool;
 
 
index 1e16a40..79cd30d 100755 (executable)
@@ -35,7 +35,7 @@
 
 #define   IS_C2_COMPOSITE 0
 //============================================================
-inline GeomAbs_CurveType Geom2dInt_CurveTool::TheType(const IntCurveCurve& C) {
+inline GeomAbs_CurveType Geom2dInt_CurveTool::GetType(const IntCurveCurve& C) {
   return(C.GetType());
 }
 //============================================================
@@ -85,6 +85,25 @@ inline void Geom2dInt_CurveTool::D2 (const IntCurveCurve& C,
   
  C.D2(U,P,T,N);
 }
+
+//============================================================
+inline void Geom2dInt_CurveTool::D3 (const IntCurveCurve& C,
+                                    const Standard_Real U,
+                                    gp_Pnt2d& P,
+                                    gp_Vec2d& T,
+                                    gp_Vec2d& N,
+                                    gp_Vec2d& V) {
+  
+ C.D3(U,P,T,N,V);
+}
+//============================================================
+inline gp_Vec2d Geom2dInt_CurveTool::DN(const Adaptor2d_Curve2d& C,
+             const Standard_Real U,
+             const Standard_Integer N)
+  {
+  return C.DN(U,N);  
+  }
+
 //============================================================
 inline Standard_Real Geom2dInt_CurveTool::FirstParameter (const IntCurveCurve& C) {
   return(C.FirstParameter());
@@ -134,4 +153,7 @@ inline Standard_Integer Geom2dInt_CurveTool::NbIntervals(const IntCurveCurve& C)
 }
 //============================================================
 
-
+ inline Standard_Integer Geom2dInt_CurveTool::Degree(const IntCurveCurve& C)
+{
+  return C.Degree();
+}
index 2036593..31869ee 100755 (executable)
@@ -188,7 +188,19 @@ is
       --  Warning : It used only for C2 aproximation
    returns Boolean
    is private; 
-
+   
+   RotateTrihedron(me;
+      Tangent    : out  Vec  from  gp; 
+      Normal     : out  Vec  from  gp; 
+      BiNormal   : out  Vec  from  gp;
+      NewTangent : in   Vec  from  gp)
+   ---Purpose: revolves the trihedron (which is determined
+   -- of given "Tangent", "Normal" and "BiNormal" vectors)
+   -- to coincide "Tangent" and "NewTangent" axes.
+   returns Boolean from Standard
+   is private; 
+   
+   
 fields 
    P         :  Pnt  from  gp;
    mySngl    :  HArray1OfReal from TColStd;  
index 964f004..8212e12 100755 (executable)
 #include <TCollection_CompareOfReal.hxx>
 #include <TColgp_SequenceOfPnt2d.hxx>
 
-#define NullTol 1.e-10
-#define MaxSingular 1.e-5
+static const Standard_Real NullTol = 1.e-10;
+static const Standard_Real MaxSingular = 1.e-5;
+
+static const Standard_Integer maxDerivOrder = 3;
+
 
 //=======================================================================
 //function : FDeriv
@@ -60,6 +63,30 @@ static gp_Vec DDeriv(const gp_Vec& F, const gp_Vec& DF, const gp_Vec& D2F)
 }
 
 //=======================================================================
+//function : CosAngle
+//purpose  : Return a cosine between vectors theV1 and theV2.
+//=======================================================================
+static Standard_Real CosAngle(const gp_Vec& theV1, const gp_Vec& theV2)
+  {
+  const Standard_Real aTol = gp::Resolution();
+
+  const Standard_Real m1 = theV1.Magnitude(), m2 = theV2.Magnitude();
+  if((m1 <= aTol) || (m2 <= aTol)) //Vectors are codirectional
+    return 1.0;
+
+  Standard_Real aCAng = theV1.Dot(theV2)/(m1*m2);
+
+  if(aCAng > 1.0)
+    aCAng = 1.0;
+
+  if(aCAng < -1.0)
+    aCAng = -1.0;
+
+
+  return aCAng;
+  }
+
+//=======================================================================
 //function : GeomFill_Frenet
 //purpose  : 
 //=======================================================================
@@ -315,35 +342,204 @@ Handle(GeomFill_TrihedronLaw) GeomFill_Frenet::Copy() const
 }
 
 //=======================================================================
+//function :  RotateTrihedron
+//purpose  :  This function revolves the trihedron (which is determined of
+//            given "Tangent", "Normal" and "BiNormal" vectors)
+//            to coincide "Tangent" and "NewTangent" axes.
+//=======================================================================
+Standard_Boolean 
+          GeomFill_Frenet::RotateTrihedron( gp_Vec& Tangent,
+                                            gp_Vec& Normal,
+                                            gp_Vec& BiNormal,
+                                            const gp_Vec& NewTangent) const
+  {
+  const Standard_Real anInfCOS = cos(Precision::Angular()); //0.99999995
+  const Standard_Real aTol = gp::Resolution();
+
+  gp_Vec anAxis = Tangent.Crossed(NewTangent);
+  const Standard_Real NT = anAxis.Magnitude();
+  if(NT <= aTol)
+    //No rotation required
+    return Standard_True;
+  else
+    anAxis /= NT;   //Normalization
+
+  const Standard_Real aPx = anAxis.X(), aPy = anAxis.Y(), aPz = anAxis.Z();
+  const Standard_Real aCAng = CosAngle(Tangent,NewTangent); //cosine
+
+  const Standard_Real anAddCAng = 1.0 - aCAng;
+  const Standard_Real aSAng = sqrt(1.0 - aCAng*aCAng);  //sine
+
+//According to rotate direction, sine of rotation angle might be 
+//positive or negative.
+//We can research to choose necessary sign. But I think, it is more
+//effectively, to rotate "Tangent" vector in both direction. After that
+//we can choose necessary rotation direction in depend of results.
+
+  const gp_Vec aV11(anAddCAng*aPx*aPx+aCAng,      anAddCAng*aPx*aPy-aPz*aSAng,  anAddCAng*aPx*aPz+aPy*aSAng);
+  const gp_Vec aV12(anAddCAng*aPx*aPx+aCAng,      anAddCAng*aPx*aPy+aPz*aSAng,  anAddCAng*aPx*aPz-aPy*aSAng);
+  const gp_Vec aV21(anAddCAng*aPx*aPy+aPz*aSAng,  anAddCAng*aPy*aPy+aCAng,      anAddCAng*aPy*aPz-aPx*aSAng);
+  const gp_Vec aV22(anAddCAng*aPx*aPy-aPz*aSAng,  anAddCAng*aPy*aPy+aCAng,      anAddCAng*aPy*aPz+aPx*aSAng);
+  const gp_Vec aV31(anAddCAng*aPx*aPz-aPy*aSAng,  anAddCAng*aPy*aPz+aPx*aSAng,  anAddCAng*aPz*aPz+aCAng);
+  const gp_Vec aV32(anAddCAng*aPx*aPz+aPy*aSAng,  anAddCAng*aPy*aPz-aPx*aSAng,  anAddCAng*aPz*aPz+aCAng);
+
+  gp_Vec aT1(Tangent.Dot(aV11), Tangent.Dot(aV21), Tangent.Dot(aV31));
+  gp_Vec aT2(Tangent.Dot(aV12), Tangent.Dot(aV22), Tangent.Dot(aV32));
+
+  if(CosAngle(aT1,NewTangent) >= CosAngle(aT2,NewTangent))
+    {
+    Tangent = aT1;
+    Normal = gp_Vec(Normal.Dot(aV11), Normal.Dot(aV21), Normal.Dot(aV31));
+    BiNormal = gp_Vec(BiNormal.Dot(aV11), BiNormal.Dot(aV21), BiNormal.Dot(aV31));
+    }
+  else
+    {
+    Tangent = aT2;
+    Normal = gp_Vec(Normal.Dot(aV12), Normal.Dot(aV22), Normal.Dot(aV32));
+    BiNormal = gp_Vec(BiNormal.Dot(aV12), BiNormal.Dot(aV22), BiNormal.Dot(aV32));
+    }
+
+  return CosAngle(Tangent,NewTangent) >= anInfCOS;
+  }
+
+
+//=======================================================================
 //function : D0
 //purpose  : 
 //=======================================================================
 
- Standard_Boolean GeomFill_Frenet::D0(const Standard_Real Param,
+ Standard_Boolean GeomFill_Frenet::D0(const Standard_Real theParam,
                                       gp_Vec& Tangent,
                                       gp_Vec& Normal,
                                       gp_Vec& BiNormal)
 {
+  const Standard_Real aTol = gp::Resolution();
+
   Standard_Real norm;
   Standard_Integer Index;
   Standard_Real Delta = 0.;
-  if(IsSingular(Param, Index)) 
-    if (SingularD0(Param, Index, Tangent, Normal, BiNormal, Delta))
+  if(IsSingular(theParam, Index)) 
+    if (SingularD0(theParam, Index, Tangent, Normal, BiNormal, Delta))
       return Standard_True;
 
-  Standard_Real theParam = Param + Delta;
-  myTrimmed->D2(theParam, P, Tangent, BiNormal);
-  Tangent.Normalize();
-  BiNormal = Tangent.Crossed(BiNormal);
-  norm = BiNormal.Magnitude();
-  if (norm <= gp::Resolution()) {
-    gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent);
-    BiNormal.SetXYZ(Axe.YDirection().XYZ());    
-  }
-  else BiNormal.Normalize();
+  Standard_Real aParam = theParam + Delta;
+  myTrimmed->D2(aParam, P, Tangent, BiNormal);
 
-  Normal = BiNormal;
-  Normal.Cross(Tangent);
+  const Standard_Real DivisionFactor = 1.e-3;
+  const Standard_Real anUinfium   = myTrimmed->FirstParameter();
+  const Standard_Real anUsupremum = myTrimmed->LastParameter();
+  const Standard_Real aDelta = (anUsupremum - anUinfium)*DivisionFactor;
+  Standard_Real Ndu = Tangent.Magnitude();
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////
+  if(Ndu <= aTol)
+    {
+    gp_Vec aTn;
+//Derivative is approximated by Taylor-series
+    
+    Standard_Integer anIndex = 1; //Derivative order
+    Standard_Boolean isDeriveFound = Standard_False;
+    
+    do
+      {
+      aTn =  myTrimmed->DN(theParam,++anIndex);
+      Ndu = aTn.Magnitude();
+      isDeriveFound = Ndu > aTol;
+      }
+    while(!isDeriveFound && anIndex < maxDerivOrder);
+
+    if(isDeriveFound)
+      {
+      Standard_Real u;
+      
+      if(theParam-anUinfium < aDelta)
+        u = theParam+aDelta;
+      else
+        u = theParam-aDelta;
+      
+      gp_Pnt P1, P2;
+      myTrimmed->D0(Min(theParam, u),P1);
+      myTrimmed->D0(Max(theParam, u),P2);
+      
+      gp_Vec V1(P1,P2);
+      Standard_Real aDirFactor = aTn.Dot(V1);
+      
+      if(aDirFactor < 0.0)
+        aTn = -aTn;
+      }//if(IsDeriveFound)
+    else
+      {
+//Derivative is approximated by three points
+
+      gp_Pnt Ptemp(0.0,0.0,0.0); //(0,0,0)-coordinate
+      gp_Pnt P1, P2, P3;
+      Standard_Boolean IsParameterGrown;
+                
+      if(theParam-anUinfium < 2*aDelta)
+        {
+        myTrimmed->D0(theParam,P1);
+        myTrimmed->D0(theParam+aDelta,P2);
+        myTrimmed->D0(theParam+2*aDelta,P3);
+        IsParameterGrown = Standard_True;
+        }
+      else
+        {
+        myTrimmed->D0(theParam-2*aDelta,P1);
+        myTrimmed->D0(theParam-aDelta,P2);
+        myTrimmed->D0(theParam,P3);
+        IsParameterGrown = Standard_False;
+        }
+        
+      gp_Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
+      
+      if(IsParameterGrown)
+        aTn=-3*V1+4*V2-V3;
+      else
+        aTn=V1-4*V2+3*V3;
+      }//else of "if(IsDeriveFound)" condition
+      Ndu = aTn.Magnitude();
+      gp_Pnt Pt = P;
+      Standard_Real dPar = 10.0*aDelta;
+
+//Recursive calling is used for determine of trihedron for 
+//point, which is near to given.
+      if(theParam-anUinfium < dPar)
+        {
+        if(D0(aParam+dPar,Tangent,Normal,BiNormal) == Standard_False)
+          return Standard_False;
+        }
+      else
+        {
+        if(D0(aParam-dPar,Tangent,Normal,BiNormal) == Standard_False)
+          return Standard_False;
+        }
+
+      P = Pt;
+
+      if(RotateTrihedron(Tangent,Normal,BiNormal,aTn) == Standard_False)
+        {
+#ifdef DEB
+        cout << "Cannot coincide two tangents." << endl;
+#endif
+        return Standard_False;
+        }
+    }//if(Ndu <= aTol)
+  else
+    {
+    Tangent = Tangent/Ndu;
+    BiNormal = Tangent.Crossed(BiNormal);
+    norm = BiNormal.Magnitude();
+    if (norm <= gp::Resolution())
+      {
+      gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent);
+      BiNormal.SetXYZ(Axe.YDirection().XYZ());    
+      }
+    else
+      BiNormal.Normalize();
+
+    Normal = BiNormal;
+    Normal.Cross(Tangent);
+    }
 
   return Standard_True;
 }
index ceb14ac..67ec846 100755 (executable)
@@ -99,6 +99,27 @@ is
         --- Purpose : Raised if the continuity of the current interval
         --  is not C2.
      is redefined static;
+
+    D3 (me; U : Real; P : out Pnt from gp; V1, V2, V3 : out Vec from gp)
+        --- Purpose :
+        --  Returns the point P of parameter U, the first, the second 
+        --  and the third derivative.
+     raises 
+       DomainError from Standard
+        --- Purpose : Raised if the continuity of the current interval
+        --  is not C1.
+
+    is redefined static;
+        
+    DN (me; U : Real; N : Integer)   returns Vec from gp
+        --- Purpose :
+        --  The returned vector gives the value of the derivative for the 
+        --  order of derivation N.
+    raises  
+        OutOfRange from Standard
+        --- Purpose : Raised if N < 1.            
+    is redefined static;
+
  
     Resolution(me; R3d : Real) returns Real
          ---Purpose :  Returns the parametric  resolution corresponding
index 859ea9c..de1bedf 100755 (executable)
@@ -21,6 +21,7 @@
 
 
 #include <GeomFill_SnglrFunc.ixx>
+#include <Standard_NotImplemented.hxx>
 #include <Precision.hxx>
 
 GeomFill_SnglrFunc::GeomFill_SnglrFunc(const Handle(Adaptor3d_HCurve)& HC) : 
@@ -121,6 +122,45 @@ void GeomFill_SnglrFunc::SetRatio(const Standard_Real Ratio)
   V2 *= ratio;
 }
 
+void GeomFill_SnglrFunc::D3(const Standard_Real U,gp_Pnt& P,gp_Vec& V1,gp_Vec& V2,gp_Vec& V3) const
+  {
+  gp_Vec DC, D2C, D3C, D4C, D5C;
+  myHCurve->D3(U, P, DC, D2C, D3C);
+  D4C = myHCurve->DN(U, 4);
+  D5C = myHCurve->DN(U, 5);
+  P = gp_Pnt(DC.Crossed(D2C).XYZ()).ChangeCoord()*ratio;
+  V1 = DC.Crossed(D3C)*ratio;
+  V2 = (D2C.Crossed(D3C) + DC.Crossed(D4C))*ratio;
+  V3 = (DC.Crossed(D5C) + D2C.Crossed(D4C)*2)*ratio;
+  }
+
+gp_Vec GeomFill_SnglrFunc::DN(const Standard_Real U,const Standard_Integer N) const
+  {
+  Standard_RangeError_Raise_if (N < 1, "Exception: Geom2d_OffsetCurve::DN(). N<1.");
+
+  gp_Vec D1C, D2C, D3C;
+  gp_Pnt C;
+
+  switch(N)
+    {
+    case 1:
+      D1(U,C,D1C);
+      return D1C;
+    case 2:
+      D2(U,C,D1C,D2C);
+      return D2C;
+    case 3:
+      D3(U,C,D1C,D2C,D3C);
+      return D3C;
+    default:
+      Standard_NotImplemented::Raise("Exception: Derivative order is greater than 3. "
+        "Cannot compute of derivative.");
+    }
+   
+  return gp_Vec();
+
+  }
+
  Standard_Real GeomFill_SnglrFunc::Resolution(const Standard_Real R3D) const
 {
   return Precision::Parametric(R3D);
index d5fe1ec..6cb85a2 100755 (executable)
@@ -55,12 +55,12 @@ Standard_IMPORT Draw_Viewer dout;
 //=======================================================================
 
 static Standard_Integer proj (Draw_Interpretor& di, Standard_Integer n, const char** a)
-{
-  if ( n < 5)
   {
+  if ( n < 5)
+    {
     cout << " Use proj curve/surf x y z [extrema algo: g(grad)/t(tree)]" << endl;
     return 1;
-  }
+    }
 
   gp_Pnt P(Draw::Atof(a[2]),Draw::Atof(a[3]),Draw::Atof(a[4]));
 
@@ -69,70 +69,88 @@ static Standard_Integer proj (Draw_Interpretor& di, Standard_Integer n, const ch
   Handle(Geom_Curve) GC = DrawTrSurf::GetCurve(a[1]);
   Handle(Geom_Surface) GS;
   Extrema_ExtAlgo aProjAlgo = Extrema_ExtAlgo_Grad;
+
   if (n == 6 && a[5][0] == 't')
     aProjAlgo = Extrema_ExtAlgo_Tree;
 
-  if (GC.IsNull())  {
+  if (GC.IsNull())
+    {
     GS = DrawTrSurf::GetSurface(a[1]);
+    
     if (GS.IsNull())
       return 1;
+
     Standard_Real U1, U2, V1, V2;
     GS->Bounds(U1,U2,V1,V2);
 
     GeomAPI_ProjectPointOnSurf proj(P,GS,U1,U2,V1,V2,aProjAlgo);
+
     Standard_Real UU,VV;
-    for ( Standard_Integer i = 1; i <= proj.NbPoints(); i++) {
+    for ( Standard_Integer i = 1; i <= proj.NbPoints(); i++)
+      {
       gp_Pnt P1 = proj.Point(i);
-      if ( P.Distance(P1) > Precision::Confusion()) {
-       Handle(Geom_Line) L = new Geom_Line(P,gp_Vec(P,P1));
-       Handle(Geom_TrimmedCurve) CT = 
-         new Geom_TrimmedCurve(L, 0., P.Distance(P1));
+      if ( P.Distance(P1) > Precision::Confusion())
+        {
+        Handle(Geom_Line) L = new Geom_Line(P,gp_Vec(P,P1));
+        Handle(Geom_TrimmedCurve) CT = 
+          new Geom_TrimmedCurve(L, 0., P.Distance(P1));
        Sprintf(name,"%s%d","ext_",i);
-       char* temp = name; // portage WNT
-       DrawTrSurf::Set(temp, CT);
-       di << name << " ";
-      }
-      else {
+        char* temp = name; // portage WNT
+        DrawTrSurf::Set(temp, CT);
+        di << name << " ";
+        }
+      else
+        {
        Sprintf(name,"%s%d","ext_",i);
-       di << name << " ";
-       char* temp = name; // portage WNT
-       DrawTrSurf::Set(temp, P1);
-       proj.Parameters(i,UU,VV);
-       di << " Le point est sur la surface." << "\n";
-       di << " Ses parametres sont:  UU = " << UU << "\n";
-       di << "                       VV = " << VV << "\n";
+        di << name << " ";
+        char* temp = name; // portage WNT
+        DrawTrSurf::Set(temp, P1);
+        proj.Parameters(i,UU,VV);
+        di << " Le point est sur la surface." << "\n";
+        di << " Ses parametres sont:  UU = " << UU << "\n";
+        di << "                       VV = " << VV << "\n";
+        }
       }
     }
-
-  }
-  else {
+  else
+    {
     GeomAPI_ProjectPointOnCurve proj(P,GC,GC->FirstParameter(),
-                                         GC->LastParameter());
-//    Standard_Real UU;
-    for ( Standard_Integer i = 1; i <= proj.NbPoints(); i++) {
+                                                GC->LastParameter());
+
+    if(proj.NbPoints() == 0)
+      {
+      cout << "No project point was found." << endl;
+      return 0;
+      }
+
+    for ( Standard_Integer i = 1; i <= proj.NbPoints(); i++)
+      {
       gp_Pnt P1 = proj.Point(i);
       Standard_Real UU = proj.Parameter(i);
       di << " parameter " << i << " = " << UU << "\n";
-      if ( P.Distance(P1) > Precision::Confusion()) {
-       Handle(Geom_Line) L = new Geom_Line(P,gp_Vec(P,P1));
-       Handle(Geom_TrimmedCurve) CT = 
-         new Geom_TrimmedCurve(L, 0., P.Distance(P1));
+      if ( P.Distance(P1) > Precision::Confusion())
+        {
+        Handle(Geom_Line) L = new Geom_Line(P,gp_Vec(P,P1));
+        Handle(Geom_TrimmedCurve) CT = 
+                      new Geom_TrimmedCurve(L, 0., P.Distance(P1));
        Sprintf(name,"%s%d","ext_",i);
-       char* temp = name; // portage WNT
-       DrawTrSurf::Set(temp, CT);
-       di << name << " ";
-      }
-      else {
+        char* temp = name; // portage WNT
+        DrawTrSurf::Set(temp, CT);
+        di << name << " ";
+        }
+      else
+        {
        Sprintf(name,"%s%d","ext_",i);
-       char* temp = name; // portage WNT
-       DrawTrSurf::Set(temp, P1);
-       di << name << " ";
-       UU = proj.Parameter(i);
-       di << " Le point est sur la courbe." << "\n";
-       di << " Son parametre est U = " << UU << "\n";
+        char* temp = name; // portage WNT
+        DrawTrSurf::Set(temp, P1);
+        di << name << " ";
+        UU = proj.Parameter(i);
+        di << " Le point est sur la courbe." << "\n";
+        di << " Son parametre est U = " << UU << "\n";
+        }
       }
     }
-  }
+
   return 0;
 }
 
@@ -349,78 +367,133 @@ static Standard_Integer extrema(Draw_Interpretor& di, Standard_Integer n, const
   }
 
   char name[100];
-  if ( C1 && C2) {
+  if ( C1 && C2)
+    {
     GeomAPI_ExtremaCurveCurve Ex(GC1,GC2,U1f,U1l,U2f,U2l);
-    if(!Ex.Extrema().IsParallel()) {
-      for ( Standard_Integer i = 1; i <= Ex.NbExtrema(); i++) {
-       gp_Pnt P1,P2;
-       Ex.Points(i,P1,P2);
-       if (P1.Distance(P2) < 1.e-16) {
-         di << "Extrema " << i << " is point : " << P1.X() << " " << P1.Y() << " " << P1.Z() << "\n";
-         continue;
-       }
-       Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2));
-       Handle(Geom_TrimmedCurve) CT = 
-         new Geom_TrimmedCurve(L, 0., P1.Distance(P2));
-       Sprintf(name,"%s%d","ext_",i);
-       char* temp = name; // portage WNT
-       DrawTrSurf::Set(temp, CT);
-       di << name << " ";
+
+    if(!Ex.Extrema().IsParallel())
+      {
+      const Standard_Integer aNExtr = Ex.NbExtrema();
+      if(aNExtr == 0)
+        {
+        di << "No solutions!\n";
+        }
+      else
+        {
+        for ( Standard_Integer i = 1; i <= aNExtr; i++)
+          {
+          gp_Pnt P1,P2;
+          Ex.Points(i,P1,P2);
+          Standard_Real U1,V1;
+          Ex.Parameters(i,U1,V1);
+          if (P1.Distance(P2) < 1.e-16)
+            {
+            di << "Extrema " << i << " is point : " << P1.X() << " " << P1.Y() << " " << P1.Z() << "\n";
+            continue;
+            }
+
+          Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2));
+          Handle(Geom_TrimmedCurve) CT = 
+            new Geom_TrimmedCurve(L, 0., P1.Distance(P2));
+#ifdef DEB          
+          Sprintf(name,"%s%d (U=%f; V=%f)","ext_",i,U1,V1);
+#else
+          Sprintf(name,"%s%d","ext_",i);
+#endif
+          char* temp = name; // portage WNT
+          DrawTrSurf::Set(temp, CT);
+          di << name << " ";
+          }
+        }
       }
-    }
-    else {
+    else
+      {
       di << "Infinite number of extremas, distance = " << Ex.LowerDistance() << "\n";
+      }
     }
-  }
-  else if ( C1 & S2) {
+  else if ( C1 & S2)
+    {
     GeomAPI_ExtremaCurveSurface Ex(GC1,GS2,U1f,U1l,U2f,U2l,V2f,V2l);
-    for ( Standard_Integer i = 1; i <= Ex.NbExtrema(); i++) {
-      gp_Pnt P1,P2;
-      Ex.Points(i,P1,P2);
-      if (P1.Distance(P2) < 1.e-16) continue;
-      Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2));
-      Handle(Geom_TrimmedCurve) CT = 
-       new Geom_TrimmedCurve(L, 0., P1.Distance(P2));
-      Sprintf(name,"%s%d","ext_",i);
-      char* temp = name; // portage WNT
-      DrawTrSurf::Set(temp, CT);
-      di << name << " ";
+
+    const Standard_Integer aNExtr = Ex.NbExtrema();
+    if(aNExtr == 0)
+      {
+      di << "No solutions!\n";
+      }
+    else
+      {
+      for ( Standard_Integer i = 1; i <= aNExtr; i++)
+        {
+        gp_Pnt P1,P2;
+        Ex.Points(i,P1,P2);
+        if (P1.Distance(P2) < 1.e-16) continue;
+        Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2));
+        Handle(Geom_TrimmedCurve) CT = 
+          new Geom_TrimmedCurve(L, 0., P1.Distance(P2));
+        Sprintf(name,"%s%d","ext_",i);
+        char* temp = name; // portage WNT
+        DrawTrSurf::Set(temp, CT);
+        di << name << " ";
+        }
+      }
     }
-  }
-  else if ( S1 & C2) {
+  else if ( S1 & C2)
+    {
     GeomAPI_ExtremaCurveSurface Ex(GC2,GS1,U2f,U2l,U1f,U1l,V1f,V1l);
-    for ( Standard_Integer i = 1; i <= Ex.NbExtrema(); i++) {
-      gp_Pnt P1,P2;
-      Ex.Points(i,P1,P2);
-      if (P1.Distance(P2) < 1.e-16) continue;
-      Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2));
-      Handle(Geom_TrimmedCurve) CT = 
-       new Geom_TrimmedCurve(L, 0., P1.Distance(P2));
-      Sprintf(name,"%s%d","ext_",i);
-      char* temp = name; // portage WNT
-      DrawTrSurf::Set(temp, CT);
-      di << name << " ";
+
+    const Standard_Integer aNExtr = Ex.NbExtrema();
+    if(aNExtr == 0)
+      {
+      di << "No solutions!\n";
+      }
+    else
+      {
+      for ( Standard_Integer i = 1; i <= aNExtr; i++)
+        {
+        gp_Pnt P1,P2;
+        Ex.Points(i,P1,P2);
+        if (P1.Distance(P2) < 1.e-16) continue;
+        Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2));
+        Handle(Geom_TrimmedCurve) CT = 
+          new Geom_TrimmedCurve(L, 0., P1.Distance(P2));
+        Sprintf(name,"%s%d","ext_",i);
+        char* temp = name; // portage WNT
+        DrawTrSurf::Set(temp, CT);
+        di << name << " ";
+        }
+      }
     }
-  }
-  else if ( S1 & S2) {
+  else if ( S1 & S2)
+    {
     GeomAPI_ExtremaSurfaceSurface Ex(GS1,GS2,U1f,U1l,V1f,V1l,U2f,U2l,V2f,V2l);
-    for ( Standard_Integer i = 1; i <= Ex.NbExtrema(); i++) {
-      gp_Pnt P1,P2;
-      Ex.Points(i,P1,P2);
-      if (P1.Distance(P2) < 1.e-16) continue;
-      Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2));
-      Handle(Geom_TrimmedCurve) CT = 
-       new Geom_TrimmedCurve(L, 0., P1.Distance(P2));
-      Sprintf(name,"%s%d","ext_",i);
-      char* temp = name; // portage WNT
-      DrawTrSurf::Set(temp, CT);
-      di << name << " ";
+
+    const Standard_Integer aNExtr = Ex.NbExtrema();
+    if(aNExtr == 0)
+      {
+      di << "No solutions!\n";
+      }
+    else
+      {
+      for ( Standard_Integer i = 1; i <= aNExtr; i++)
+        {
+        gp_Pnt P1,P2;
+        Ex.Points(i,P1,P2);
+        if (P1.Distance(P2) < 1.e-16)
+          continue;
+
+        Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2));
+        Handle(Geom_TrimmedCurve) CT = 
+          new Geom_TrimmedCurve(L, 0., P1.Distance(P2));
+        Sprintf(name,"%s%d","ext_",i);
+        char* temp = name; // portage WNT
+        DrawTrSurf::Set(temp, CT);
+        di << name << " ";
+        }
+      }
     }
-  }
 
   return 0;
-
-}
+  }
 
 //=======================================================================
 //function : totalextcc
index fe4878b..c130e0f 100755 (executable)
@@ -1329,150 +1329,219 @@ static Standard_Integer surfpoints (Draw_Interpretor& /*di*/, Standard_Integer /
 //purpose  : 
 //=======================================================================
 static Standard_Integer intersection (Draw_Interpretor& di, Standard_Integer n, const char** a)
-{
-  if (n < 4) {
+  {
+  if (n < 4)
     return 1;
-  }
+  
   //
   Handle(Geom_Curve) GC1;
   Handle(Geom_Surface) GS1 = DrawTrSurf::GetSurface(a[2]);
-  if (GS1.IsNull()) {
+  if (GS1.IsNull())
+    {
     GC1 = DrawTrSurf::GetCurve(a[2]);
     if (GC1.IsNull())
       return 1;
-  }
+    }
+
   //
   Handle(Geom_Surface) GS2 = DrawTrSurf::GetSurface(a[3]);
-  if (GS2.IsNull()) {
+  if (GS2.IsNull())
     return 1;
-  }
+
   //
   Standard_Real tol = Precision::Confusion();
-  if (n == 5 || n == 9 || n == 13 || n == 17) tol = Draw::Atof(a[n-1]);
+  if (n == 5 || n == 9 || n == 13 || n == 17)
+    tol = Draw::Atof(a[n-1]);
+
   //
   Handle(Geom_Curve) Result;
   gp_Pnt             Point;
+  
   //
-  if (GC1.IsNull()) {
+  if (GC1.IsNull())
+    {
     GeomInt_IntSS Inters;
     //
     // Surface Surface
-    if (n <= 5) {
+    if (n <= 5)
+      {
       // General case
       Inters.Perform(GS1,GS2,tol,Standard_True);
-    }
-    else if (n == 8 || n == 9 || n == 12 || n == 13 || n == 16 || n == 17) {
+      }
+    else if (n == 8 || n == 9 || n == 12 || n == 13 || n == 16 || n == 17)
+      {
       Standard_Boolean useStart = Standard_True, useBnd = Standard_True;
       Standard_Integer ista1=0,ista2=0,ibnd1=0,ibnd2=0;
       Standard_Real UVsta[4];
       Handle(GeomAdaptor_HSurface) AS1,AS2;
+      
       //
-      if (n <= 9) {        // user starting point
+      if (n <= 9)          // user starting point
+        {
         useBnd = Standard_False;
-        ista1 = 4; ista2 = 7;
-      }
-      else if (n <= 13) {        // user bounding
+        ista1 = 4;
+        ista2 = 7;
+        }
+      else if (n <= 13)   // user bounding
+        {
         useStart = Standard_False;
         ibnd1 = 4; ibnd2 = 11;
-      }
-      else {        // both user starting point and bounding
+        }
+      else        // both user starting point and bounding
+        {
         ista1 = 4; ista2 = 7;
         ibnd1 = 8; ibnd2 = 15;
-      }
+        }
+
       if (useStart)
+        {
         for (Standard_Integer i=ista1; i <= ista2; i++)
+          {
           UVsta[i-ista1] = Draw::Atof(a[i]);
-      if (useBnd) {
+          }
+        }
+
+      if (useBnd)
+        {
         Standard_Real UVbnd[8];
         for (Standard_Integer i=ibnd1; i <= ibnd2; i++)
           UVbnd[i-ibnd1] = Draw::Atof(a[i]);
+
         AS1 = new GeomAdaptor_HSurface(GS1,UVbnd[0],UVbnd[1],UVbnd[2],UVbnd[3]);
         AS2 = new GeomAdaptor_HSurface(GS2,UVbnd[4],UVbnd[5],UVbnd[6],UVbnd[7]);
-      }
+        }
+
       //
-      if (useStart && !useBnd) {
+      if (useStart && !useBnd)
+        {
         Inters.Perform(GS1,GS2,tol,UVsta[0],UVsta[1],UVsta[2],UVsta[3]);
-      }
-      else if (!useStart && useBnd) {
+        }
+      else if (!useStart && useBnd)
+        {
         Inters.Perform(AS1,AS2,tol);
-      }
-      else {
+        }
+      else
+        {
         Inters.Perform(AS1,AS2,tol,UVsta[0],UVsta[1],UVsta[2],UVsta[3]);
-      }
-    }//else if (n == 8 || n == 9 || n == 12 || n == 13 || n == 16 || n == 17) {
-    else {
+        }
+      }//else if (n == 8 || n == 9 || n == 12 || n == 13 || n == 16 || n == 17)
+    else
+      {
       di<<"incorrect number of arguments"<<"\n";
       return 1;
-    }
+      }
+
     //
-    if (!Inters.IsDone()) {
+    if (!Inters.IsDone())
+      {
+      di<<"No intersections found!"<<"\n";
+
       return 1;
-    }
+      }
+    
     //
     char buf[1024];
     Standard_Integer i, aNbLines, aNbPoints; 
-      //
+    
+    //
     aNbLines = Inters.NbLines();
-    if (aNbLines >= 2) {
-      for (i=1; i<=aNbLines; ++i) {
-       Sprintf(buf, "%s_%d",a[1],i);
-       Result = Inters.Line(i);
-       const char* temp = buf; 
-       DrawTrSurf::Set(temp,Result);
+    if (aNbLines >= 2)
+      {
+      for (i=1; i<=aNbLines; ++i)
+        {
+        Sprintf(buf, "%s_%d",a[1],i);
+        di << buf << " ";
+        Result = Inters.Line(i);
+        const char* temp = buf;
+        DrawTrSurf::Set(temp,Result);
+        }
       }
-    }
-    else if (aNbLines == 1) {
+    else if (aNbLines == 1)
+      {
       Result = Inters.Line(1);
+      Sprintf(buf,"%s",a[1]);
+      di << buf << " ";
       DrawTrSurf::Set(a[1],Result);
     }
+
     //
     aNbPoints=Inters.NbPoints();
-    for (i=1; i<=aNbPoints; ++i) {
+    for (i=1; i<=aNbPoints; ++i)
+      {
       Point=Inters.Point(i);
       Sprintf(buf,"%s_p_%d",a[1],i);
-      const char* temp =buf;
+      di << buf << " ";
+      const char* temp = buf;
       DrawTrSurf::Set(temp, Point);
-    }
-  }// if (GC1.IsNull()) {
-  //
-  else {
+      }
+    }// if (GC1.IsNull())
+  else
+    {
     // Curve Surface
     GeomAPI_IntCS Inters(GC1,GS2);
-    //
-    if (!Inters.IsDone()) return 1;
     
+    //
+    if (!Inters.IsDone())
+      {
+      di<<"No intersections found!"<<"\n";
+      return 1;
+      }
+
     Standard_Integer nblines = Inters.NbSegments();
     Standard_Integer nbpoints = Inters.NbPoints();
-    if ( (nblines+nbpoints) >= 2) {
-      char newname[1024];
+
+    char newname[1024];
+
+    if ( (nblines+nbpoints) >= 2)
+      {
       Standard_Integer i;
       Standard_Integer Compt = 1;
-      for (i = 1; i <= nblines; i++, Compt++) {
-       Sprintf(newname,"%s_%d",a[1],Compt);
-       Result = Inters.Segment(i);
-       const char* temp = newname; // pour portage WNT
-       DrawTrSurf::Set(temp,Result);
-      }
-      for (i = 1; i <= nbpoints; i++, Compt++) {
-       Sprintf(newname,"%s_%d",a[1],i);
-       Point = Inters.Point(i);
-       const char* temp = newname; // pour portage WNT
-       DrawTrSurf::Set(temp,Point);
+
+      if(nblines >= 1)
+        cout << "   Lines: " << endl;
+
+      for (i = 1; i <= nblines; i++, Compt++)
+        {
+        Sprintf(newname,"%s_%d",a[1],Compt);
+        di << newname << " ";
+        Result = Inters.Segment(i);
+        const char* temp = newname; // pour portage WNT
+        DrawTrSurf::Set(temp,Result);
+        }
+
+      if(nbpoints >= 1)
+        cout << "   Points: " << endl;
+
+      const Standard_Integer imax = nblines+nbpoints;
+
+      for (/*i = 1*/; i <= imax; i++, Compt++)
+        {
+        Sprintf(newname,"%s_%d",a[1],i);
+        di << newname << " ";
+        Point = Inters.Point(i);
+        const char* temp = newname; // pour portage WNT
+        DrawTrSurf::Set(temp,Point);
+        }
       }
-    }
-    else if (nblines == 1) {
+    else if (nblines == 1)
+      {
       Result = Inters.Segment(1);
+      Sprintf(newname,"%s",a[1]);
+      di << newname << " ";
       DrawTrSurf::Set(a[1],Result);
-    }
-    else if (nbpoints == 1) {
+      }
+    else if (nbpoints == 1)
+      {
       Point = Inters.Point(1);
+      Sprintf(newname,"%s",a[1]);
+      di << newname << " ";
       DrawTrSurf::Set(a[1],Point);
+      }
     }
-  }
 
   dout.Flush();
   return 0;
-}
+  }
 
 //=======================================================================
 //function : CurveCommands
index f8eeeb6..83b40f1 100755 (executable)
@@ -245,7 +245,8 @@ static Standard_Integer extrema(Draw_Interpretor& di, Standard_Integer n, const
 // modified by APV (compilation error - LINUX)
 //  for ( Standard_Integer i = 1; i <= Ex.NbExtrema(); i++) {
   Standard_Integer i;
-  for ( i = 1; i <= Ex.NbExtrema(); i++) {
+  const Standard_Integer aNExtr = Ex.NbExtrema();
+  for ( i = 1; i <= aNExtr; i++) {
 // modified by APV (compilation error - LINUX)
 
     gp_Pnt2d P1,P2;
@@ -267,7 +268,7 @@ static Standard_Integer extrema(Draw_Interpretor& di, Standard_Integer n, const
     }
   }
   if (i==1)
-    di << "No decisions ";
+    di << "No solutions!\n";
 
   return 0;
 }
index 918d4fb..003c396 100755 (executable)
@@ -245,4 +245,8 @@ is
     NbSamples(myclass; C: Address from Standard)
     returns Integer from Standard;
 
+    Degree(myclass; C: Address from Standard)
+      returns Integer from Standard;
+      ---C++: inline
+
 end CurveTool;
index b5ed07f..2065e82 100755 (executable)
@@ -313,3 +313,15 @@ inline Handle(Geom2d_BSplineCurve)
 inline Standard_Real
 HLRBRep_CurveTool::EpsX(const Standard_Address C)
 { return(1e-10); }
+
+
+//=======================================================================
+//function : Degree
+//purpose  : 
+//=======================================================================
+
+inline Standard_Integer
+     HLRBRep_CurveTool::Degree (const Standard_Address C)
+{
+  return(((HLRBRep_Curve *)C)->Degree());
+}
index 08003db..9273f42 100755 (executable)
@@ -39,7 +39,7 @@ void IntCurve_IntCurveCurveGen::Perform(const TheCurve& C,
   IntRes2d_Domain D1;
   Standard_Real TolDomain = Tol;
   if(Tol<TolConf) TolDomain = TolConf;
-  GeomAbs_CurveType typ = TheCurveTool::TheType(C);
+  GeomAbs_CurveType typ = TheCurveTool::GetType(C);
   switch(typ) {        
   case GeomAbs_Ellipse:        
   case GeomAbs_Circle:
@@ -95,7 +95,7 @@ void IntCurve_IntCurveCurveGen::Perform(const TheCurve& C,
                                        const IntRes2d_Domain& D,
                                        const Standard_Real TolConf,
                                        const Standard_Real Tol) { 
-  GeomAbs_CurveType typ = TheCurveTool::TheType(C);
+  GeomAbs_CurveType typ = TheCurveTool::GetType(C);
   switch(typ) {        
   case GeomAbs_Ellipse:        
   case GeomAbs_Circle:
@@ -125,7 +125,7 @@ IntRes2d_Domain IntCurve_IntCurveCurveGen::ComputeDomain(const TheCurve& C1,
                                                         const Standard_Real TolDomain) const { 
   IntRes2d_Domain D1;
 
-  GeomAbs_CurveType typ = TheCurveTool::TheType(C1);
+  GeomAbs_CurveType typ = TheCurveTool::GetType(C1);
   switch(typ) { 
     
   case GeomAbs_Ellipse: 
@@ -238,8 +238,8 @@ void IntCurve_IntCurveCurveGen::InternalPerform (const TheCurve& C1,
                                                 const Standard_Real Tol,
                                                 const Standard_Boolean Composite) {
 
-  GeomAbs_CurveType typ1 = TheCurveTool::TheType(C1);
-  GeomAbs_CurveType typ2 = TheCurveTool::TheType(C2);
+  GeomAbs_CurveType typ1 = TheCurveTool::GetType(C1);
+  GeomAbs_CurveType typ2 = TheCurveTool::GetType(C2);
 
 
   switch (typ1) {
index 9b02927..caf7c28 100755 (executable)
@@ -324,7 +324,7 @@ void IntCurve_UserIntConicCurveGen::InternalPerform (const gp_Lin2d& Lin1,
                                                     const Standard_Real Tol,
                                                     const Standard_Boolean Composite) {
 
-  GeomAbs_CurveType typ2 = ThePCurveTool::TheType(C2);
+  GeomAbs_CurveType typ2 = ThePCurveTool::GetType(C2);
   
   switch (typ2) {
     
@@ -419,7 +419,7 @@ void IntCurve_UserIntConicCurveGen::InternalPerform (const gp_Circ2d& Circ1,
                                                     const Standard_Real Tol,
                                                     const Standard_Boolean Composite) { 
 
-  GeomAbs_CurveType typ2 = ThePCurveTool::TheType(C2);
+  GeomAbs_CurveType typ2 = ThePCurveTool::GetType(C2);
   
   switch (typ2) {
     
@@ -513,7 +513,7 @@ void IntCurve_UserIntConicCurveGen::InternalPerform (const gp_Elips2d& Elips1,
                                                     const Standard_Real Tol,
                                                     const Standard_Boolean Composite) { 
   
-  GeomAbs_CurveType typ2 = ThePCurveTool::TheType(C2);
+  GeomAbs_CurveType typ2 = ThePCurveTool::GetType(C2);
     
   switch (typ2) {
     
@@ -609,7 +609,7 @@ void IntCurve_UserIntConicCurveGen::InternalPerform (const gp_Parab2d& Parab1,
                                                     const Standard_Real Tol,
                                                     const Standard_Boolean Composite) { 
 
-  GeomAbs_CurveType typ2 = ThePCurveTool::TheType(C2);
+  GeomAbs_CurveType typ2 = ThePCurveTool::GetType(C2);
   
   switch (typ2) {
     
@@ -704,7 +704,7 @@ void IntCurve_UserIntConicCurveGen::InternalPerform (const gp_Hypr2d& Hyper1,
                                                     const Standard_Real Tol,
                                                     const Standard_Boolean Composite) { 
 
-  GeomAbs_CurveType typ2 = ThePCurveTool::TheType(C2);
+  GeomAbs_CurveType typ2 = ThePCurveTool::GetType(C2);
   
   switch (typ2) {
     
index 69b1e41..9af5890 100755 (executable)
@@ -139,21 +139,21 @@ is
 
 fields
 
-    myCurve   : Curve;    -- the Curve on which thw calculus are done
-    u         : Real;     -- the current value of the parameter
-    level     : Integer;  -- the order of derivation
-    cn        : Real;     -- the order of continuity of the Curve
-    linTol    : Real;     -- the tolerance for null Vector
+    myCurve      : Curve;    -- the Curve on which thw calculus are done
+    myU          : Real;     -- the current value of the parameter
+    myDerOrder   : Integer;  -- the order of derivation
+    myCN         : Real;     -- the order of continuity of the Curve
+    myLinTol     : Real;     -- the tolerance for null Vector
     
-    pnt       : Pnt;      -- the current point value
-    d         : Vec[3];   -- the current first, second and third derivative
-                          -- value
-    tangent   : Dir;      -- the tangent value
-    curvature : Real;     -- the curvature value
+    myPnt        : Pnt;      -- the current point value
+    myDerivArr   : Vec[3];   -- the current first, second and third derivative
+                             -- value
+    myTangent    : Dir;      -- the tangent value
+    myCurvature  : Real;     -- the curvature value
 
-    tangentStatus   : Status from LProp; 
+    myTangentStatus   : Status from LProp; 
                           -- the status of the tangent direction
-    significantFirstDerivativeOrder : Integer;
+    mySignificantFirstDerivativeOrder : Integer;
                           -- the order of the first non null derivative
                           -- 
 end CLProps;
index 14a38ad..b4819d3 100755 (executable)
 #include <LProp_NotDefined.hxx>
 #include <Standard_OutOfRange.hxx>
 
+static const Standard_Real MinStep   = 1.0e-7;
 
 
-LProp_CLProps::LProp_CLProps (const Curve& C, 
+
+LProp_CLProps::LProp_CLProps (const Curve& C,
                               const Standard_Real U,
-                             const Standard_Integer N, 
+                              const Standard_Integer N, 
                               const Standard_Real Resolution) 
-     : myCurve(C), 
-       level(N),
-       cn(4), // cn(Tool::Continuity(C)),   RLE
-       linTol(Resolution), 
-       tangentStatus (LProp_Undecided)
+      : myCurve(C), myDerOrder(N), myCN(4), 
+        myLinTol(Resolution), myTangentStatus (LProp_Undecided)
 {
-
   Standard_OutOfRange_Raise_if (N < 0 || N > 3,
-                                "LProp_CLProps::LProp_CLProps()");
-  
+                          "LProp_CLProps::LProp_CLProps()");
+
   SetParameter(U);
 }
 
-LProp_CLProps::LProp_CLProps (const Curve& C, 
-                             const Standard_Integer N, 
-                              const Standard_Real Resolution) 
-     : myCurve(C), 
-       u(RealLast()),
-       level(N),
-       cn(4), // (Tool::Continuity(C)), RLE
-       linTol(Resolution), 
-       tangentStatus (LProp_Undecided)
-
+LProp_CLProps::LProp_CLProps (const Curve& C, const Standard_Integer N, 
+                                    const Standard_Real Resolution) 
+        : myCurve(C), myU(RealLast()), myDerOrder(N), myCN(4), 
+            myLinTol(Resolution), myTangentStatus (LProp_Undecided)
 {
-
-  Standard_OutOfRange_Raise_if (N < 0 || N > 3, "");  
+  Standard_OutOfRange_Raise_if (N < 0 || N > 3, 
+                          "LProp_CLProps::LProp_CLProps()");  
 }
 
 LProp_CLProps::LProp_CLProps (const Standard_Integer N, 
                               const Standard_Real Resolution) 
-     : u(RealLast()),
-       level(N),
-       cn(0),
-       linTol(Resolution),
-       tangentStatus (LProp_Undecided)
+        : myU(RealLast()), myDerOrder(N), myCN(0), myLinTol(Resolution),
+                            myTangentStatus (LProp_Undecided)
 {
-
   Standard_OutOfRange_Raise_if (N < 0 || N > 3, "");  
 }
 
 void LProp_CLProps::SetParameter(const Standard_Real U)
-{ 
-  u = U;
-  switch (level) {
+{
+  myU = U;
+  switch (myDerOrder)
+  {
   case 0:
-    Tool::Value(myCurve, u, pnt);
+    Tool::Value(myCurve, myU, myPnt);
     break;
   case 1:
-    Tool::D1(myCurve, u, pnt, d[0]);
+    Tool::D1(myCurve, myU, myPnt, myDerivArr[0]);
     break;
   case 2:
-    Tool::D2(myCurve, u, pnt, d[0], d[1]);
+    Tool::D2(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1]);
     break;
   case 3:
-    Tool::D3(myCurve, u, pnt, d[0], d[1], d[2]);
+    Tool::D3(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1], myDerivArr[2]);
     break;
   }
-  tangentStatus = LProp_Undecided;
+
+  myTangentStatus = LProp_Undecided;
 }
 
-void LProp_CLProps::SetCurve(const Curve& C) {
-       myCurve = C ; 
-       cn = 4; // Tool::Continuity(C); RLE
+void LProp_CLProps::SetCurve(const Curve& C)
+{
+  myCurve = C ; 
+  myCN = 4; // Tool::Continuity(C); RLE
 }
 
 const Pnt& LProp_CLProps::Value () const
-{ 
-  return pnt;
+{
+  return myPnt;
 }
 
 const Vec& LProp_CLProps::D1 ()
 {
-  if (level < 1) {
-    level = 1;
-    Tool::D1(myCurve, u, pnt, d[0]);
+  if (myDerOrder < 1)
+  {
+    myDerOrder = 1;
+    Tool::D1(myCurve, myU, myPnt, myDerivArr[0]);
   }
-  return d[0];
+
+  return myDerivArr[0];
 }
 
 const Vec& LProp_CLProps::D2 ()
 {
-  if (level < 2) {
-    level = 2;
-    Tool::D2(myCurve, u, pnt, d[0], d[1]);
+  if (myDerOrder < 2)
+  {
+    myDerOrder = 2;
+    Tool::D2(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1]);
   }
-  return d[1];
+
+  return myDerivArr[1];
 }
 
 const Vec& LProp_CLProps::D3 ()
 {
-  if (level < 3) {
-    level = 3;
-    Tool::D3(myCurve, u, pnt, d[0], d[1], d[2]);
+  if (myDerOrder < 3)
+  {
+    myDerOrder = 3;
+    Tool::D3(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1], myDerivArr[2]);
   }
-  return d[2];
+
+  return myDerivArr[2];
 }
 
 Standard_Boolean LProp_CLProps::IsTangentDefined ()
 {
-
-  if (tangentStatus == LProp_Undefined) {
+  if (myTangentStatus == LProp_Undefined)
     return Standard_False;
-  }
-  else if (tangentStatus >= LProp_Defined) {
+  else if (myTangentStatus >= LProp_Defined)
     return Standard_True;
-  }
 
   // tangentStatus == Lprop_Undecided 
   // we have to calculate the first non null derivative
-  Standard_Real Tol = linTol * linTol;
+  const Standard_Real Tol = myLinTol * myLinTol;
+  
   Vec V;
+  
   Standard_Integer Order = 0;
-  while (Order < 4) {
-    Order++;
-    if(cn >= Order) {
-      switch(Order) {
-      case 1 :
-       V = D1();
-       break;
-      case 2 :
-       V = D2();
-       break;
-      case 3 :
-       V = D3();
-       break;
-      };
-      if(V.SquareMagnitude() > Tol) {
-       significantFirstDerivativeOrder = Order;
-       tangentStatus = LProp_Defined;
-       return Standard_True;
-      }
-    }
-    else {
-      tangentStatus = LProp_Undefined;
+  while (Order++ < 4)
+  {
+    if(myCN >= Order)
+    {
+      switch(Order)
+      {
+      case 1:
+        V = D1();
+        break;
+      case 2:
+        V = D2();
+        break;
+      case 3:
+        V = D3();
+        break;
+      }//switch(Order)
+
+      if(V.SquareMagnitude() > Tol)
+      {
+        mySignificantFirstDerivativeOrder = Order;
+        myTangentStatus = LProp_Defined;
+        return Standard_True;
+      }//if(V.SquareMagnitude() > Tol)
+    }//if(cn >= Order)
+    else
+    {
+      myTangentStatus = LProp_Undefined;
       return Standard_False;
-    }
-  }
+    }// else of "if(cn >= Order)" condition
+  }//while (Order < 4)
+
   return Standard_False;
 }
 
-
 void  LProp_CLProps::Tangent (Dir& D)
 {
-
-  if(!IsTangentDefined()) { LProp_NotDefined::Raise(); }
-  D = Dir(d[significantFirstDerivativeOrder - 1]);
+  if(!IsTangentDefined())
+    LProp_NotDefined::Raise();
+  
+  if(mySignificantFirstDerivativeOrder == 1)
+    D = Dir(myDerivArr[0]);
+  else if (mySignificantFirstDerivativeOrder > 1)
+  {
+    const Standard_Real DivisionFactor = 1.e-3;
+    const Standard_Real anUsupremum = Tool::LastParameter(myCurve), 
+                        anUinfium = Tool::FirstParameter(myCurve);
+                        
+    Standard_Real du;
+    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst()))
+      du = 0.0;
+    else
+      du = anUsupremum-anUinfium;
+    
+    const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
+
+    Vec V = myDerivArr[mySignificantFirstDerivativeOrder - 1];
+    
+    Standard_Real u;
+          
+    if(myU-anUinfium < aDelta)
+      u = myU+aDelta;
+    else
+      u = myU-aDelta;
+    
+    Pnt P1, P2;
+    Tool::Value(myCurve, Min(myU, u),P1);
+    Tool::Value(myCurve, Max(myU, u),P2);
+    
+    Vec V1(P1,P2);
+    Standard_Real aDirFactor = V.Dot(V1);
+    
+    if(aDirFactor < 0.0)
+      V = -V;
+      
+    D = Dir(V);
+  }//else if (mySignificantFirstDerivativeOrder > 1)
 }
 
 Standard_Real LProp_CLProps::Curvature ()
 {
   Standard_Boolean isDefined = IsTangentDefined();
   LProp_NotDefined_Raise_if(!isDefined,
-                           "LProp_CLProps::CurvatureNotDefined()");
+          "LProp_CLProps::CurvatureNotDefined()");
 
   // if the first derivative is null the curvature is infinite.
-  if(significantFirstDerivativeOrder > 1) return RealLast();
+  if(mySignificantFirstDerivativeOrder > 1)
+    return RealLast();
 
-  Standard_Real Tol = linTol * linTol;
-  Standard_Real DD1 = d[0].SquareMagnitude();
-  Standard_Real DD2 = d[1].SquareMagnitude();
+  Standard_Real Tol = myLinTol * myLinTol;
+  Standard_Real DD1 = myDerivArr[0].SquareMagnitude();
+  Standard_Real DD2 = myDerivArr[1].SquareMagnitude();
+  
   // if the second derivative is null the curvature is null.
-  if (DD2 <= Tol) { 
-    curvature = 0.0;
+  if (DD2 <= Tol)
+  {
+    myCurvature = 0.0;
   }
-  else {
-    Standard_Real N = d[0].CrossSquareMagnitude(d[1]);
+  else
+  {
+    Standard_Real N = myDerivArr[0].CrossSquareMagnitude(myDerivArr[1]);
     // if d[0] and d[1] are colinear the curvature is null.
     Standard_Real t = N/(DD1*DD2);
-    if (t<=Tol) {
-      curvature = 0.0;
+    if (t<=Tol)
+    {
+      myCurvature = 0.0;
     }
-    else {
-      curvature = sqrt(N) / (DD1*sqrt(DD1));
+    else
+    {
+      myCurvature = sqrt(N) / (DD1*sqrt(DD1));
     }
   }
-  return curvature;
-} 
 
+  return myCurvature;
+}
 
 void  LProp_CLProps::Normal (Dir& D)
 {
   Standard_Real c = Curvature();
-  if(c==RealLast() || Abs(c) <= linTol) { LProp_NotDefined::Raise(); }
+  if(c==RealLast() || Abs(c) <= myLinTol)
+  {
+    LProp_NotDefined::Raise("LProp_CLProps::Normal(...):"
+        "Curvature is null or infinity"); 
+  }
 
   // we used here the following vector relation 
   // a ^ (b ^ c) = b(ac) - c(ab)
   // Norm = d[0] ^ (d[1] ^ d[0])
-
-  Vec Norm = d[1] * (d[0] * d[0]) - d[0] * (d[0] * d[1]);
+  
+  Vec Norm = myDerivArr[1] * (myDerivArr[0] * myDerivArr[0]) - myDerivArr[0] * (myDerivArr[0] * myDerivArr[1]);
   D = Dir(Norm);
 }
 
-
 void  LProp_CLProps::CentreOfCurvature (Pnt& P)
 {
-  if(Abs(Curvature()) <= linTol) { LProp_NotDefined::Raise(); }
+  if(Abs(Curvature()) <= myLinTol)
+  {
+    LProp_NotDefined::Raise(); 
+  }
 
   // we used here the following vector relation 
   // a ^ (b ^ c) = b(ac) - c(ab)
   // Norm = d[0] ^ (d[1] ^ d[0])
 
-  Vec Norm = d[1] * (d[0] * d[0]) - d[0] * (d[0] * d[1]);
+  Vec Norm = myDerivArr[1] * (myDerivArr[0] * myDerivArr[0]) - myDerivArr[0] * (myDerivArr[0] * myDerivArr[1]);
   Norm.Normalize();
-  Norm.Divide(curvature);
-  P= pnt.Translated(Norm);
+  Norm.Divide(myCurvature);
+  P= myPnt.Translated(Norm);
 }
index 3e30651..8ae9190 100755 (executable)
@@ -197,35 +197,35 @@ is
 
 fields
 
-    surf   : Surface;
-    u      : Real;
-    v      : Real;
-    level  : Integer;
-    cn     : Integer;
-    linTol : Real;
-
-    pnt  : Pnt from gp;
-    d1U  : Vec from gp;
-    d1V  : Vec from gp;
-    d2U  : Vec from gp;
-    d2V  : Vec from gp;
-    dUV  : Vec from gp;
-
-    normal     : Dir from gp;
-    minCurv    : Real;
-    maxCurv    : Real;
-    dirMinCurv : Dir from gp;
-    dirMaxCurv : Dir from gp;
-    meanCurv   : Real;
-    gausCurv   : Real;
-
-    significantFirstUDerivativeOrder : Integer;
-    significantFirstVDerivativeOrder : Integer;
-
-    uTangentStatus  : Status from LProp;
-    vTangentStatus  : Status from LProp;
-    normalStatus    : Status from LProp;
-    curvatureStatus : Status from LProp;
+    mySurf     : Surface;
+    myU        : Real;
+    myV        : Real;
+    myDerOrder  : Integer;
+    myCN       : Integer;
+    myLinTol   : Real;
+
+    myPnt  : Pnt from gp;
+    myD1u  : Vec from gp;
+    myD1v  : Vec from gp;
+    myD2u  : Vec from gp;
+    myD2v  : Vec from gp;
+    myDuv  : Vec from gp;
+
+    myNormal     : Dir from gp;
+    myMinCurv    : Real;
+    myMaxCurv    : Real;
+    myDirMinCurv : Dir from gp;
+    myDirMaxCurv : Dir from gp;
+    myMeanCurv   : Real;
+    myGausCurv   : Real;
+
+    mySignificantFirstDerivativeOrderU : Integer;
+    mySignificantFirstDerivativeOrderV : Integer;
+
+    myUTangentStatus  : Status from LProp;
+    myVTangentStatus  : Status from LProp;
+    myNormalStatus    : Status from LProp;
+    myCurvatureStatus : Status from LProp;
 
 end SLProps;
 
index 66604b6..608522e 100755 (executable)
 #include <TColgp_Array2OfVec.hxx>
 #include <math_DirectPolynomialRoots.hxx>
 
+static const Standard_Real MinStep   = 1.0e-7;
 
 static Standard_Boolean IsTangentDefined (LProp_SLProps& SProp,
-                                 const Standard_Integer  cn,
-                                 const Standard_Real     linTol,
-                                 const Standard_Integer  Derivative, 
-                                 Standard_Integer&       Order,
-                                 LProp_Status&  Status) 
+                                          const Standard_Integer cn,
+                                          const Standard_Real linTol,
+                                          const Standard_Integer  Derivative, 
+                                          Standard_Integer&       Order,
+                                          LProp_Status&  Status) 
 {
   Standard_Real Tol = linTol * linTol;
   gp_Vec V[2];
   Order = 0;
-  while (Order < 3) {
+
+  while (Order < 3)
+  {
     Order++;
-    if(cn >= Order) {
-      switch(Order) {
-      case 1 :
-       V[0] = SProp.D1U();
-       V[1] = SProp.D1V();
-       break;
-      case 2 :
-       V[0] = SProp.D2U();
-       V[1] = SProp.D2V();
-       break;
-      };
-      if(V[Derivative].SquareMagnitude() > Tol) {
-       Status = LProp_Defined;
-       return Standard_True;
+    if(cn >= Order)
+    {
+      switch(Order)
+      {
+      case 1:
+        V[0] = SProp.D1U();
+        V[1] = SProp.D1V();
+        break;
+      case 2:
+        V[0] = SProp.D2U();
+        V[1] = SProp.D2V();
+        break;
+      }//switch(Order)
+
+      if(V[Derivative].SquareMagnitude() > Tol)
+      {
+        Status = LProp_Defined;
+        return Standard_True;
       }
-    }
-    else {
+    }//if(cn >= Order)
+    else
+    {
       Status = LProp_Undefined;
       return Standard_False;
     }
   }
+
   return Standard_False;
 }
 
 LProp_SLProps::LProp_SLProps (const Surface& S, 
-                             const Standard_Real     U,  
-                              const Standard_Real     V, 
-                              const Standard_Integer  N, 
-                              const Standard_Real     Resolution) 
-     : surf(S),
-       level(N),
-       cn(4), // (Tool::Continuity(S)),
-       linTol(Resolution)
+                              const Standard_Real U,  
+                              const Standard_Real V, 
+                              const Standard_Integer N, 
+                              const Standard_Real Resolution) 
+        : mySurf(S),myDerOrder(N), myCN(4), // (Tool::Continuity(S)),
+            myLinTol(Resolution)
 {
-
   Standard_OutOfRange_Raise_if(N < 0 || N > 2, 
-                               "LProp_SLProps::LProp_SLProps()");
+                        "LProp_SLProps::LProp_SLProps()");
 
   SetParameters(U, V);
 }
@@ -83,292 +89,328 @@ LProp_SLProps::LProp_SLProps (const Surface& S,
 LProp_SLProps::LProp_SLProps (const Surface& S, 
                               const Standard_Integer  N, 
                               const Standard_Real     Resolution) 
-     : surf(S),
-       u(RealLast()), v(RealLast()),
-       level(N),
-       cn(4), // (Tool::Continuity(S)),
-       linTol(Resolution),
-       uTangentStatus (LProp_Undecided),
-       vTangentStatus (LProp_Undecided),
-       normalStatus   (LProp_Undecided),
-       curvatureStatus(LProp_Undecided)
+        : mySurf(S), myU(RealLast()), myV(RealLast()), myDerOrder(N),
+            myCN(4), // (Tool::Continuity(S))
+            myLinTol(Resolution), 
+            myUTangentStatus (LProp_Undecided),
+            myVTangentStatus (LProp_Undecided),
+            myNormalStatus   (LProp_Undecided),
+            myCurvatureStatus(LProp_Undecided)
 {
   Standard_OutOfRange_Raise_if(N < 0 || N > 2, 
-                               "LProp_SLProps::LProp_SLProps()");
+                      "LProp_SLProps::LProp_SLProps()");
 }
 
 LProp_SLProps::LProp_SLProps (const Standard_Integer  N, 
-                              const Standard_Real     Resolution) 
-       :u(RealLast()), v(RealLast()),
-       level(N),
-       cn(0),
-       linTol(Resolution),
-       uTangentStatus (LProp_Undecided),
-       vTangentStatus (LProp_Undecided),
-       normalStatus   (LProp_Undecided),
-       curvatureStatus(LProp_Undecided)
+                              const Standard_Real Resolution) 
+        : myU(RealLast()), myV(RealLast()), myDerOrder(N), myCN(0),
+        myLinTol(Resolution),
+        myUTangentStatus (LProp_Undecided),
+        myVTangentStatus (LProp_Undecided),
+        myNormalStatus   (LProp_Undecided),
+        myCurvatureStatus(LProp_Undecided)
 {
   Standard_OutOfRange_Raise_if(N < 0 || N > 2, 
-                               "LProp_SLProps::LProp_SLProps() bad level");
+                "LProp_SLProps::LProp_SLProps() bad level");
 }
 
-void LProp_SLProps::SetSurface (const Surface& S ) {
-
-     surf = S;
-     cn = 4; // =Tool::Continuity(S);
+void LProp_SLProps::SetSurface (const Surface& S ) 
+{
+  mySurf = S;
+  myCN = 4; // =Tool::Continuity(S);
 }
 
-void LProp_SLProps::SetParameters (const Standard_Real U, const Standard_Real V)
+void LProp_SLProps::SetParameters (const Standard_Real U, 
+                                   const Standard_Real V)
 {
-  u = U;
-  v = V;
-  switch (level) {
+  myU = U;
+  myV = V;
+  switch (myDerOrder)
+  {
   case 0:
-    Tool::Value(surf, u, v, pnt);
+    Tool::Value(mySurf, myU, myV, myPnt);
     break;
   case 1:
-    Tool::D1(surf, u, v, pnt, d1U, d1V);
+    Tool::D1(mySurf, myU, myV, myPnt, myD1u, myD1v);
     break;
   case 2:
-    Tool::D2(surf, u, v, pnt, d1U, d1V, d2U, d2V, dUV);
+    Tool::D2(mySurf, myU, myV, myPnt, myD1u, myD1v, myD2u, myD2v, myDuv);
     break;
-  };
-  uTangentStatus  = LProp_Undecided;
-  vTangentStatus  = LProp_Undecided;
-  normalStatus    = LProp_Undecided;
-  curvatureStatus = LProp_Undecided;
+  }
+
+  myUTangentStatus  = LProp_Undecided;
+  myVTangentStatus  = LProp_Undecided;
+  myNormalStatus    = LProp_Undecided;
+  myCurvatureStatus = LProp_Undecided;
 }
 
 const gp_Pnt& LProp_SLProps::Value() const
 {
-  return pnt;
+  return myPnt;
 }
 
 const gp_Vec& LProp_SLProps::D1U()
 {
-  if (level < 1) {
-    level =1;
-    Tool::D1(surf,u,v,pnt,d1U,d1V);
+  if (myDerOrder < 1)
+  {
+    myDerOrder =1;
+    Tool::D1(mySurf,myU,myV,myPnt,myD1u,myD1v);
   }
-  return d1U;
+
+  return myD1u;
 }
 
 const gp_Vec& LProp_SLProps::D1V()
 {
-  if (level < 1) {
-    level =1;
-    Tool::D1(surf,u,v,pnt,d1U,d1V);
+  if (myDerOrder < 1)
+  {
+    myDerOrder =1;
+    Tool::D1(mySurf,myU,myV,myPnt,myD1u,myD1v);
   }
-  return d1V;
+
+  return myD1v;
 }
 
 const gp_Vec& LProp_SLProps::D2U()
 {
-  if (level < 2) {
-    level =2;
-    Tool::D2(surf,u,v,pnt,d1U,d1V,d2U,d2V,dUV);
+  if (myDerOrder < 2) 
+  {
+    myDerOrder =2;
+    Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv);
   }
-  return d2U;
+
+  return myD2u;
 }
 
 const gp_Vec& LProp_SLProps::D2V()
 {
-  if (level < 2) {
-    level =2;
-    Tool::D2(surf,u,v,pnt,d1U,d1V,d2U,d2V,dUV);
+  if (myDerOrder < 2) 
+  {
+    myDerOrder =2;
+    Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv);
   }
-  return d2V;
+
+  return myD2v;
 }
+
 const gp_Vec& LProp_SLProps::DUV()
 {
-  if (level < 2) {
-    level =2;
-    Tool::D2(surf,u,v,pnt,d1U,d1V,d2U,d2V,dUV);
+  if (myDerOrder < 2) 
+  {
+    myDerOrder =2;
+    Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv);
   }
-  return dUV;
+
+  return myDuv;
 }
 
 Standard_Boolean LProp_SLProps::IsTangentUDefined ()
 {
-  if (uTangentStatus == LProp_Undefined) { 
+  if (myUTangentStatus == LProp_Undefined)
     return Standard_False;
-  }
-  else if (uTangentStatus >= LProp_Defined) { 
+  else if (myUTangentStatus >= LProp_Defined)
     return Standard_True; 
-  }
+
   // uTangentStatus == Lprop_Undecided 
   // we have to calculate the first non null U derivative
-  return IsTangentDefined(*this, cn, linTol, 0, 
-                          significantFirstUDerivativeOrder, uTangentStatus);
+  return IsTangentDefined(*this, myCN, myLinTol, 0, 
+          mySignificantFirstDerivativeOrderU, myUTangentStatus);
 }
 
-void LProp_SLProps::TangentU (gp_Dir& D) {
-
-  if(!IsTangentUDefined()) { LProp_NotDefined::Raise(); }
-  if(significantFirstUDerivativeOrder == 1) {
-    D = gp_Dir(d1U);
-  }
-  else {
-    D = gp_Dir(d2U);
+void LProp_SLProps::TangentU (gp_Dir& D) 
+{
+  if(!IsTangentUDefined())
+    LProp_NotDefined::Raise();
+
+  if(mySignificantFirstDerivativeOrderU == 1)
+    D = gp_Dir(myD1u);
+  else
+  {
+    const Standard_Real DivisionFactor = 1.e-3;
+    Standard_Real anUsupremum, anUinfium;
+    Standard_Real anVsupremum, anVinfium;
+    Tool::Bounds(mySurf,anUinfium,anVinfium,anUsupremum,anVsupremum);
+    Standard_Real du;
+    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst()))
+      du = 0.0;
+    else
+      du = anUsupremum-anUinfium;
+    
+    const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep);
+
+    gp_Vec V = myD2u;
+    
+    Standard_Real u;
+          
+    if(myU-anUinfium < aDeltaU)
+      u = myU+aDeltaU;
+    else
+      u = myU-aDeltaU;
+    
+    gp_Pnt P1, P2;
+    Tool::Value(mySurf, Min(myU, u),myV,P1);
+    Tool::Value(mySurf, Max(myU, u),myV,P2);
+    
+    gp_Vec V1(P1,P2);
+    Standard_Real aDirFactor = V.Dot(V1);
+    
+    if(aDirFactor < 0.0)
+      V = -V;
+    
+    D = gp_Dir(V);
   }
 }
 
 Standard_Boolean LProp_SLProps::IsTangentVDefined ()
 {
-
-  if (vTangentStatus == LProp_Undefined) { 
+  if (myVTangentStatus == LProp_Undefined)
     return Standard_False;
-  }
-  else if (vTangentStatus >= LProp_Defined) { 
+  else if (myVTangentStatus >= LProp_Defined)
     return Standard_True; 
-  }
+
   // vTangentStatus == Lprop_Undecided 
   // we have to calculate the first non null V derivative
-  return IsTangentDefined(*this, cn, linTol, 1, 
-                          significantFirstVDerivativeOrder, vTangentStatus);
+  return IsTangentDefined(*this, myCN, myLinTol, 1, 
+          mySignificantFirstDerivativeOrderV, myVTangentStatus);
 }
 
-void LProp_SLProps::TangentV (gp_Dir& D) {
-
-  if(!IsTangentVDefined()) { LProp_NotDefined::Raise(); }
-  if(significantFirstVDerivativeOrder == 1) {
-    D = gp_Dir(d1V);
-  }
-  else {
-    D = gp_Dir(d2V);
+void LProp_SLProps::TangentV (gp_Dir& D)
+{
+  if(!IsTangentVDefined())
+    LProp_NotDefined::Raise();
+  
+  if(mySignificantFirstDerivativeOrderV == 1)
+    D = gp_Dir(myD1v);
+  else
+  {
+    const Standard_Real DivisionFactor = 1.e-3;
+    Standard_Real anUsupremum, anUinfium;
+    Standard_Real anVsupremum, anVinfium;
+    Tool::Bounds(mySurf,anUinfium,anVinfium,anUsupremum,anVsupremum);
+    
+    Standard_Real dv;
+    if((anVsupremum >= RealLast()) || (anVinfium <= RealFirst()))
+      dv = 0.0;
+    else
+      dv = anVsupremum-anVinfium;
+    
+    const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep);
+
+    gp_Vec V = myD2v;
+    
+    Standard_Real v;
+          
+    if(myV-anVinfium < aDeltaV)
+      v = myV+aDeltaV;
+    else
+      v = myV-aDeltaV;
+    
+    gp_Pnt P1, P2;
+    Tool::Value(mySurf, myU, Min(myV, v),P1);
+    Tool::Value(mySurf, myU, Max(myV, v),P2);
+    
+    gp_Vec V1(P1,P2);
+    Standard_Real aDirFactor = V.Dot(V1);
+    
+    if(aDirFactor < 0.0)
+      V = -V;
+    
+    D = gp_Dir(V);
   }
 }
 
 Standard_Boolean LProp_SLProps::IsNormalDefined()
 {
-
-  if (normalStatus == LProp_Undefined) {
+  if (myNormalStatus == LProp_Undefined)
     return Standard_False;
-  }
-  else if (normalStatus >= LProp_Defined) {
+  else if (myNormalStatus >= LProp_Defined)
     return Standard_True;
-  }
+
   // status = UnDecided 
-  
   // first try the standard computation of the normal.
   CSLib_DerivativeStatus Status;
-  CSLib::Normal(d1U, d1V, linTol, Status, normal);
-  if (Status  == CSLib_Done ) { 
-    normalStatus = LProp_Computed;
+  CSLib::Normal(myD1u, myD1v, myLinTol, Status, myNormal);
+  if (Status  == CSLib_Done ) 
+  {
+    myNormalStatus = LProp_Computed;
     return Standard_True;
   }
-  // else solve the degenerated case only if continuity >= 2
-/*  if (cn >= 2) {
-    if(level < 2) this->D2U();
-    Standard_Boolean Done;
-    CSLib_NormalStatus Stat;
-    CSLib::Normal(d1U, d1V, d2U, d2V, dUV, linTol, Done, Stat, normal);
-    if (Done) {
-      normalStatus = LProp_Computed;
-      return Standard_True;
-    }
-  }*/
-  /*
-  else {
-    Standard_Integer MaxOrder=3;
-    CSLib_NormalStatus Stat;
-    gp_Dir thenormal;
-    TColgp_Array2OfVec DerNUV(0,MaxOrder,0,MaxOrder);
-    TColgp_Array2OfVec DerSurf(0,MaxOrder+1,0,MaxOrder+1);
-    Standard_Integer i,j,OrderU,OrderV;
-    Standard_Real Umin,Umax,Vmin,Vmax;
-    Tool::Bounds(surf, Umin, Vmin, Umax, Vmax);
-
-    // Calcul des derivees
-    for(i=1;i<=MaxOrder+1;i++){
-      DerSurf.SetValue(i,0, Tool::DN(surf,u,v,i,0));
-    }
-  
-    for(i=0;i<=MaxOrder+1;i++)
-      for(j=1;j<=MaxOrder+1;j++){
-       DerSurf.SetValue(i,j, Tool::DN(surf,u,v,i,j));
-      }
 
-    for(i=0;i<=MaxOrder;i++)
-      for(j=0;j<=MaxOrder;j++){
-       DerNUV.SetValue(i,j,CSLib::DNNUV(i,j,DerSurf));
-      }
-
-    CSLib::Normal(MaxOrder,DerNUV, 1.e-9, u,v,
-                 Umin,Umax,Vmin,Vmax,
-                 Stat, thenormal, OrderU, OrderV);
-    normal.SetXYZ(thenormal.XYZ());
-    if (Stat == CSLib_Defined) {
-      normalStatus = LProp_Computed;
-      return Standard_True;
-    }
-  }
-  */
+  // else solve the degenerated case only if continuity >= 2
 
-  normalStatus = LProp_Undefined;
+  myNormalStatus = LProp_Undefined;
   return Standard_False;
 }
 
-const gp_Dir& LProp_SLProps::Normal () {
+const gp_Dir& LProp_SLProps::Normal ()
+{
+  if(!IsNormalDefined())
+  {
+    LProp_NotDefined::Raise();
+  }
 
-  if(!IsNormalDefined()) { LProp_NotDefined::Raise(); }
-  return normal;
+  return myNormal;
 }
 
+
 Standard_Boolean LProp_SLProps::IsCurvatureDefined ()
 {
-  
-  if (curvatureStatus == LProp_Undefined) {
+  if (myCurvatureStatus == LProp_Undefined)
     return Standard_False;
-  }
-  else if (curvatureStatus >= LProp_Defined) {
+  else if (myCurvatureStatus >= LProp_Defined)
     return Standard_True;
-  }
-  if(cn < 2) {
-    curvatureStatus = LProp_Undefined;
+  
+  if(myCN < 2)
+  {
+    myCurvatureStatus = LProp_Undefined;
     return Standard_False;
   }
+
   // status = UnDecided 
-  if (!IsNormalDefined()) {
-    curvatureStatus = LProp_Undefined;
+  if (!IsNormalDefined())
+  {
+    myCurvatureStatus = LProp_Undefined;
     return Standard_False;
   }
+
   // pour eviter un plantage dans le cas du caro pointu
   // en fait on doit pouvoir calculer les courbure
   // avoir
-  if(!IsTangentUDefined() || !IsTangentVDefined()) {
-    curvatureStatus = LProp_Undefined;
+  if(!IsTangentUDefined() || !IsTangentVDefined())
+  {
+    myCurvatureStatus = LProp_Undefined;
     return Standard_False;
   }
 
   // here we compute the curvature features of the surface
 
-  gp_Vec Norm (normal);
+  gp_Vec Norm (myNormal);
 
-  Standard_Real E = d1U.SquareMagnitude();
-  Standard_Real F = d1U.Dot(d1V);
-  Standard_Real G = d1V.SquareMagnitude();
+  Standard_Real E = myD1u.SquareMagnitude();
+  Standard_Real F = myD1u.Dot(myD1v);
+  Standard_Real G = myD1v.SquareMagnitude();
 
-  if(level < 2) this->D2U();
+  if(myDerOrder < 2)
+    this->D2U();
 
-  Standard_Real L = Norm.Dot(d2U);
-  Standard_Real M = Norm.Dot(dUV);
-  Standard_Real N = Norm.Dot(d2V);
+  Standard_Real L = Norm.Dot(myD2u);
+  Standard_Real M = Norm.Dot(myDuv);
+  Standard_Real N = Norm.Dot(myD2v);
 
   Standard_Real A = E * M - F * L;
   Standard_Real B = E * N - G * L;
   Standard_Real C = F * N - G * M;
 
   Standard_Real MaxABC = Max(Max(Abs(A),Abs(B)),Abs(C));
-  if (MaxABC < RealEpsilon()) {          // ombilic
-    minCurv = N / G;
-    maxCurv = minCurv;
-    dirMinCurv = gp_Dir (d1U);
-    dirMaxCurv = gp_Dir (d1U.Crossed(Norm));
-    meanCurv = minCurv;            // (Cmin + Cmax) / 2.
-    gausCurv = minCurv * minCurv;  // (Cmin * Cmax)
-    curvatureStatus = LProp_Computed;
+  if (MaxABC < RealEpsilon())    // ombilic
+  {
+    myMinCurv = N / G;
+    myMaxCurv = myMinCurv;
+    myDirMinCurv = gp_Dir (myD1u);
+    myDirMaxCurv = gp_Dir (myD1u.Crossed(Norm));
+    myMeanCurv = myMinCurv;            // (Cmin + Cmax) / 2.
+    myGausCurv = myMinCurv * myMinCurv;  // (Cmin * Cmax)
+    myCurvatureStatus = LProp_Computed;
     return Standard_True;
   }
 
@@ -377,105 +419,124 @@ Standard_Boolean LProp_SLProps::IsCurvatureDefined ()
   C = C / MaxABC;
   Standard_Real Curv1, Curv2, Root1, Root2;
   gp_Vec VectCurv1, VectCurv2;
-  if (Abs(A) > RealEpsilon()) {
+
+  if (Abs(A) > RealEpsilon())
+  {
     math_DirectPolynomialRoots Root (A, B, C);
-    if(Root.NbSolutions() != 2) {
-      curvatureStatus = LProp_Undefined;
+    if(Root.NbSolutions() != 2)
+    {
+      myCurvatureStatus = LProp_Undefined;
       return Standard_False;
     }
-    else {
-       Root1 = Root.Value(1);
+    else
+    {
+      Root1 = Root.Value(1);
       Root2 = Root.Value(2);
       Curv1 = ((L * Root1 + 2. * M) * Root1 + N) /
              ((E * Root1 + 2. * F) * Root1 + G);
       Curv2 = ((L * Root2 + 2. * M) * Root2 + N) /
              ((E * Root2 + 2. * F) * Root2 + G);
-      VectCurv1 = Root1 * d1U + d1V;
-      VectCurv2 = Root2 * d1U + d1V;
+      VectCurv1 = Root1 * myD1u + myD1v;
+      VectCurv2 = Root2 * myD1u + myD1v;
     }
   }
-  else if (Abs(C) > RealEpsilon()) {
+  else if (Abs(C) > RealEpsilon())
+  {
     math_DirectPolynomialRoots Root(C, B, A);
-    if((Root.NbSolutions() != 2)) {
-      curvatureStatus = LProp_Undefined;
+
+    if((Root.NbSolutions() != 2))
+    {
+      myCurvatureStatus = LProp_Undefined;
       return Standard_False;
     }
-    else {
+    else
+    {
       Root1 = Root.Value(1);
       Root2 = Root.Value(2);
       Curv1 = ((N * Root1 + 2. * M) * Root1 + L) /
              ((G * Root1 + 2. * F) * Root1 + E);
       Curv2 = ((N * Root2 + 2. * M) * Root2 + L) /
              ((G * Root2 + 2. * F) * Root2 + E);
-      VectCurv1 = d1U + Root1 * d1V;
-      VectCurv2 = d1U + Root2 * d1V;
+      VectCurv1 = myD1u + Root1 * myD1v;
+      VectCurv2 = myD1u + Root2 * myD1v;
     }
   }
-  else {
+  else
+  {
     Curv1 = L / E;
     Curv2 = N / G;
-    VectCurv1 = d1U;
-    VectCurv2 = d1V;
+    VectCurv1 = myD1u;
+    VectCurv2 = myD1v;
   }
-  if (Curv1 < Curv2) {
-    minCurv = Curv1;
-    maxCurv = Curv2;
-    dirMinCurv = gp_Dir (VectCurv1);
-    dirMaxCurv = gp_Dir (VectCurv2);
+
+  if (Curv1 < Curv2)
+  {
+    myMinCurv = Curv1;
+    myMaxCurv = Curv2;
+    myDirMinCurv = gp_Dir (VectCurv1);
+    myDirMaxCurv = gp_Dir (VectCurv2);
   }
-  else {
-    minCurv = Curv2;
-    maxCurv = Curv1;
-    dirMinCurv = gp_Dir (VectCurv2);
-    dirMaxCurv = gp_Dir (VectCurv1);
+  else
+  {
+    myMinCurv = Curv2;
+    myMaxCurv = Curv1;
+    myDirMinCurv = gp_Dir (VectCurv2);
+    myDirMaxCurv = gp_Dir (VectCurv1);
   }
-  meanCurv = ((N * E) - (2. * M * F) + (L * G)) // voir Farin p.282
+
+  myMeanCurv = ((N * E) - (2. * M * F) + (L * G)) // voir Farin p.282
                / (2. * ((E * G) - (F * F)));
-  gausCurv = ((L * N) - (M * M)) 
+  myGausCurv = ((L * N) - (M * M)) 
                / ((E * G) - (F * F));
-  curvatureStatus = LProp_Computed;
+  myCurvatureStatus = LProp_Computed;
   return Standard_True;
-} 
-
+}
 
 Standard_Boolean LProp_SLProps::IsUmbilic ()
 {
-  if(!IsCurvatureDefined()) { LProp_NotDefined::Raise(); }
-  return Abs(maxCurv - minCurv) < Abs(Epsilon(maxCurv));
+  if(!IsCurvatureDefined())
+    LProp_NotDefined::Raise();
+  
+  return Abs(myMaxCurv - myMinCurv) < Abs(Epsilon(myMaxCurv));
 }
 
 Standard_Real LProp_SLProps::MaxCurvature ()
 {
-  if(!IsCurvatureDefined()) { LProp_NotDefined::Raise(); }
-  return maxCurv;
+  if(!IsCurvatureDefined())
+    LProp_NotDefined::Raise();
+
+  return myMaxCurv;
 }
 
 Standard_Real LProp_SLProps::MinCurvature ()
 {
-  if(!IsCurvatureDefined()) { LProp_NotDefined::Raise(); }
-  return minCurv;
+  if(!IsCurvatureDefined())
+    LProp_NotDefined::Raise();
+  
+  return myMinCurv;
 }
 
 void LProp_SLProps::CurvatureDirections(gp_Dir& Max, gp_Dir& Min)
 {
+  if(!IsCurvatureDefined())
+    LProp_NotDefined::Raise();
 
-  if(!IsCurvatureDefined()) { LProp_NotDefined::Raise(); }
-  Max = dirMaxCurv;
-  Min = dirMinCurv;
+  Max = myDirMaxCurv;
+  Min = myDirMinCurv;
 }
 
-Standard_Real LProp_SLProps::MeanCurvature () {
-
-  if(!IsCurvatureDefined()) { LProp_NotDefined::Raise(); }
-  return meanCurv;
+Standard_Real LProp_SLProps::MeanCurvature ()
+{
+  if(!IsCurvatureDefined())
+    LProp_NotDefined::Raise();
 
+  return myMeanCurv;
 }
 
-Standard_Real LProp_SLProps::GaussianCurvature () {
+Standard_Real LProp_SLProps::GaussianCurvature ()
+{
+  if(!IsCurvatureDefined())
+    LProp_NotDefined::Raise();
 
-  if(!IsCurvatureDefined()) { LProp_NotDefined::Raise(); }
-  return gausCurv;
+  return myGausCurv;
 }
-
-
-
index 6bebd43..291d592 100755 (executable)
@@ -151,7 +151,7 @@ is
      raises VectorWithNullMagnitude
       
      is static;
-
+     
   IsParallel (me; Other : Vec; AngularTolerance : Real)   returns Boolean
         --- Purpose :
         --  Returns True if Angle(<me>, Other) <= AngularTolerance or
@@ -164,7 +164,7 @@ is
      raises VectorWithNullMagnitude
      
      is static;
-
+     
   Angle (me; Other : Vec)  returns Real
         --- Purpose :
         --  Computes the angular value between <me> and <Other>
index d77d101..6d2dd32 100755 (executable)
@@ -136,7 +136,7 @@ is
      raises VectorWithNullMagnitude
       
      is static;
-
+     
   IsParallel (me; Other : Vec2d; AngularTolerance : Real)
      returns Boolean
         ---C++: inline
@@ -212,6 +212,11 @@ is
         --- Purpose : Computes the scalar product
         ---C++: alias operator *
 
+  GetNormal(me)  returns Vec2d is static;
+        ---C++: inline 
+        -- Purpose : Returns a vector {Y(), -X()} which
+        --      is normal to given.
+
   Multiply (me : in out; Scalar : Real)           is static;
         ---C++: inline
         ---C++: alias operator *=
index fa6e3e4..8585c6f 100755 (executable)
@@ -244,3 +244,8 @@ inline gp_Vec2d operator* (const Standard_Real Scalar,
                           const gp_Vec2d& V)
 { return V.Multiplied(Scalar); }
 
+inline gp_Vec2d gp_Vec2d::GetNormal() const
+  {
+  return gp_Vec2d(this->Y(), (-1)*this->X());
+  }
+
diff --git a/tests/bugs/modalg_5/bug23706_1 b/tests/bugs/modalg_5/bug23706_1
new file mode 100755 (executable)
index 0000000..cc98cbe
--- /dev/null
@@ -0,0 +1,32 @@
+puts "============"
+puts "OCC23706"
+puts "============"
+puts ""
+#########################################################################
+# Cannot project point on curve
+#########################################################################
+
+bsplinecurve r2 4 3 1 5 2 1 3 5 0 8 0 1 2 8 2 1 4 8 3 1 4 8 3 1 6 8 4 1 10 8 10 1
+mkedge spine r2
+wire spine spine
+circle profile 0 8 0 1 0 1 0.2
+mkedge profile profile
+wire profile profile
+mkplane profile profile
+pipe p spine profile
+explode p f
+mksurface ss1 p_1
+mksurface ss2 p_2
+mksurface ss3 p_3
+offset o1 ss1 0.1
+offset o2 ss2 0.1
+offset o3 ss3 0.1
+
+mkface res o2
+set info [sprops res]
+regexp {Mass +: +([-0-9.+eE]+)} $info full sq
+set sq_check 9.00819
+
+if { [expr 1.*abs($sq_check - $sq)/$sq_check] > 0.01 } {
+   puts "Error : The square of result shape is $sq"
+}
diff --git a/tests/bugs/modalg_5/bug23706_10 b/tests/bugs/modalg_5/bug23706_10
new file mode 100755 (executable)
index 0000000..491b9f0
--- /dev/null
@@ -0,0 +1,17 @@
+puts "============"
+puts "OCC23706"
+puts "============"
+puts ""
+#########################################################################
+# Cannot project point on curve
+#########################################################################
+
+bsplinecurve r3 2 6 1 3 2 1 3 1 4 1 5 1 6 3 2 5 3 1 3 7 3 1 4 8 3 1 4 8 3 1 4 8 3 1 5 9 3 1 9 7 3 1
+bsplinecurve r4 2 6 2 3 2.5 1 3 1 3.5 1 4 1 4.5 3 -1 2 3 1 1 11 3 1 3 9 3 1 3 9 3 1 3 9 3 1 5 7 3 1 7 4 3 1
+set info [extrema r3 r4]
+
+if { [regexp "Extrema 3 is point : 4 8 3" $info] != 1 || [regexp "Extrema 8 is point : 4 8 3" $info] != 1 } {
+    puts "Error : Point of extrema is wrong"
+} else {
+    puts "OK: Point of extrema is valid"
+}
diff --git a/tests/bugs/modalg_5/bug23706_11 b/tests/bugs/modalg_5/bug23706_11
new file mode 100755 (executable)
index 0000000..207023b
--- /dev/null
@@ -0,0 +1,17 @@
+puts "============"
+puts "OCC23706"
+puts "============"
+puts ""
+#########################################################################
+# Cannot project point on curve
+#########################################################################
+
+bsplinecurve r1 2 5 1 3 2 1 3 1 4 1 5 3 2 5 3 1 3 7 3 1 4 8 3 1 4 8 3 1 5 9 3 1 9 7 3 1
+bsplinecurve r2 2 5 2 3 2.5 1 3 1 3.5 1 4 3 -1 2 3 1 1 11 3 1 3 9 3 1 3 9 3 1 3 9 3 1 5 7 3 1 7 4 3 1
+set info [extrema r1 r2]
+
+if { [llength $info] != 3 } {
+    puts "Error : Extrema is wrong"
+} else {
+    puts "OK: Extrema is valid"
+}
diff --git a/tests/bugs/modalg_5/bug23706_12 b/tests/bugs/modalg_5/bug23706_12
new file mode 100755 (executable)
index 0000000..743e4e0
--- /dev/null
@@ -0,0 +1,18 @@
+puts "============"
+puts "OCC23706"
+puts "============"
+puts ""
+#########################################################################
+# Cannot project point on curve
+#########################################################################
+
+bsplinecurve r9 2 6 1 3 2 1 3 1 4 1 5 1 6 3 4 -3 3 1 6 8 3 1 10 11 3 1 10 11 3 1 10 11 3 1 14 14 3 1 5 8 3 1
+bsplinecurve r10 2 6 2 3 2.5 1 3 1 3.5 1 4 1 4.5 3 5 20 3 1 8 15 3 1 12 18 3 1 12 18 3 1 12 18 3 1 16 21 3 1 7 12 3 1
+
+set info [extrema r9 r10]
+
+if { [llength $info] != 6 } {
+    puts "Error : Extrema is wrong"
+} else {
+    puts "OK: Extrema is valid"
+}
diff --git a/tests/bugs/modalg_5/bug23706_13 b/tests/bugs/modalg_5/bug23706_13
new file mode 100755 (executable)
index 0000000..6a1f221
--- /dev/null
@@ -0,0 +1,26 @@
+puts "============"
+puts "OCC23706"
+puts "============"
+puts ""
+#########################################################################
+# Cannot project point on curve
+#########################################################################
+
+2dbsplinecurve b3 2 6 1 3 2 1 3 1 4 1 5 1 6 3 2 5 1 3 7 1 4 8 1 4 8 1 4 8 1 5 9 1 9 7 1
+2dbsplinecurve b4 2 6 2 3 2.5 1 3 1 3.5 1 4 1 4.5 3 -1 2 1 1 11 1 3 9 1 3 9 1 3 9 1 5 7 1 7 4 1
+
+set info [2dextrema b3 b4]
+set status 0
+for { set i 1 } { $i <= 15 } { incr i 1 } {
+    regexp "dist $i: +(\[-0-9.+eE\]+)" $info full pp
+    if { $pp != 1.4142135623730951 } {
+       puts "Error : Extrema is wrong on dist $i"
+       set status 1
+    }
+}
+
+if { $status != 0 } {
+    puts "Error : Extrema is wrong"
+} else {
+    puts "OK: Extrema is valid"
+}
diff --git a/tests/bugs/modalg_5/bug23706_14 b/tests/bugs/modalg_5/bug23706_14
new file mode 100755 (executable)
index 0000000..f15381d
--- /dev/null
@@ -0,0 +1,49 @@
+puts "============"
+puts "OCC23706"
+puts "============"
+puts ""
+#########################################################################
+# Cannot project point on curve
+#########################################################################
+
+2dbsplinecurve b9 2 8 1 2 2 1 3 1 4 1 5 1 6 1 7 1 8 2 4 -3 1 6 8 1 10 11 1 10 11 1 10 11 1 14 14 1 5 8 1
+2dbsplinecurve b10 2 8 2 2 2.5 1 3 1 3.5 1 4 1 4.5 1 5 1 5.5 2 5 20 1 8 15 1 12 18 1 12 18 1 12 18 1 16 21 1 7 12 1
+set info [2dextrema b9 b10]
+
+set status 0
+for { set i 1 } { $i <= 9 } { incr i } {
+    regexp "dist $i: +(\[-0-9.+eE\]+)" $info full pp1
+    if { $pp1 != 7.2801098892805181 } {
+       puts "Error : Extrema is wrong on dist $i"
+       set status 1
+    }
+}
+
+for { set j 11 } { $j <= 19 } { incr j 1 } {
+    regexp "dist $j: +(\[-0-9.+eE\]+)" $info full pp2
+    if { $pp2 != 7.2801098892805181 } {
+       puts "Error : Extrema is wrong on dist $j"
+       set status 1
+    }
+}
+
+regexp {dist 10: +([-0-9.+eE]+)} $info full pp3
+regexp {dist 20: +([-0-9.+eE]+)} $info full pp4
+regexp {dist 21: +([-0-9.+eE]+)} $info full pp5
+set pp_c 4.316921907096102
+set pp_ch 4.3169219070961038
+
+if { $pp3 != $pp_c } {
+    puts "Error : Extrema is wrong on dist 10"
+    set status 1
+}
+if { $pp4 != $pp_ch || $pp5 != $pp_ch} {
+    puts "Error : Extrema is wrong on dist 20 or 21"
+    set status 1
+}
+
+if { $status != 0 } {
+    puts "Error : Extrema is wrong"
+} else {
+    puts "OK: Extrema is valid"
+}
diff --git a/tests/bugs/modalg_5/bug23706_15 b/tests/bugs/modalg_5/bug23706_15
new file mode 100755 (executable)
index 0000000..c5a7f06
--- /dev/null
@@ -0,0 +1,34 @@
+puts "============"
+puts "OCC23706"
+puts "============"
+puts ""
+#########################################################################
+# Cannot project point on curve
+#########################################################################
+
+2dbsplinecurve b7 2 5 1 3 2 1 3 1 4 1 5 3 4 -3 1 6 8 1 10 11 1 10 11 1 14 14 1 5 8 1
+2dbsplinecurve b8 2 5 2 3 2.5 1 3 1 3.5 1 4 3 5 20 1 8 15 1 12 18 1 12 18 1 16 21 1 7 12 1
+set info [2dextrema b7 b8]
+
+set status 0
+for { set i 2 } { $i <= 5 } { incr i } {
+    regexp "dist $i: +(\[-0-9.+eE\]+)" $info full pp1
+    if { $pp1 !=4.3624023150195192 } {
+       puts "Error : Extrema is wrong on dist $i"
+       set status 1
+    }
+}
+
+regexp {dist 1: +([-0-9.+eE]+)} $info full pp2
+set pp_ch 4.3624023150195184
+
+if { $pp2 != $pp_ch } {
+    puts "Error : Extrema is wrong on dist 1"
+    set status 1
+}
+
+if { $status != 0 } {
+    puts "Error : Extrema is wrong"
+} else {
+    puts "OK: Extrema is valid"
+}
diff --git a/tests/bugs/modalg_5/bug23706_16 b/tests/bugs/modalg_5/bug23706_16
new file mode 100644 (file)
index 0000000..9fd030f
--- /dev/null
@@ -0,0 +1,22 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.999999652077201
+set y 5.0000000062915735
+set z 5.00002142991819367
+set pp_ch 0.9991079538920743
+
+restore [locate_data_file bug23706_c.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp
+if { $pp != $pp_ch } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_17 b/tests/bugs/modalg_5/bug23706_17
new file mode 100644 (file)
index 0000000..8928dff
--- /dev/null
@@ -0,0 +1,22 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.99999991301930024
+set y 5.00000000157289337
+set z 5.00000535747954842
+set pp_ch 0.99955486819730277
+
+restore [locate_data_file bug23706_c.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp
+if { $pp != $pp_ch } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_19 b/tests/bugs/modalg_5/bug23706_19
new file mode 100644 (file)
index 0000000..90a6521
--- /dev/null
@@ -0,0 +1,22 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.99999999837571056
+set y 5.0000000000293724
+set z 5.0000001000463034
+set pp_ch 0.99993927567416474
+
+restore [locate_data_file bug23706_c.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp
+if { $pp != $pp_ch } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_2 b/tests/bugs/modalg_5/bug23706_2
new file mode 100755 (executable)
index 0000000..10cc9a1
--- /dev/null
@@ -0,0 +1,26 @@
+puts "============"
+puts "OCC23706"
+puts "============"
+puts ""
+#########################################################################
+# Cannot project point on curve
+#########################################################################
+
+bsplinecurve r2 4 3 1 5 2 1 3 5 0 8 0 1 2 8 2 1 4 8 3 1 4 8 3 1 6 8 4 1 10 8 10 1
+mkedge spine r2
+wire spine spine
+mksweep spine
+addsweep spine -R
+buildsweep spine -R
+explode spine f
+mksurface ss spine_1
+offset o1 ss 2
+
+mkface res o1
+set info [sprops res]
+regexp {Mass +: +([-0-9.+eE]+)} $info full sq
+set sq_check 254.476
+
+if { [expr 1.*abs($sq_check - $sq)/$sq_check] > 0.01 } {
+   puts "Error : The square of result shape is $sq"
+}
diff --git a/tests/bugs/modalg_5/bug23706_20 b/tests/bugs/modalg_5/bug23706_20
new file mode 100644 (file)
index 0000000..a40965a
--- /dev/null
@@ -0,0 +1,22 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.9999965207720098
+set y 5.0000000629157348
+set z 5.0002142991819367
+set pp_ch 0.99715423329884789
+
+restore [locate_data_file bug23706_c2.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp
+if { $pp != $pp_ch } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_21 b/tests/bugs/modalg_5/bug23706_21
new file mode 100644 (file)
index 0000000..7f2cf8b
--- /dev/null
@@ -0,0 +1,22 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.999999652077201
+set y 5.0000000062915735
+set z 5.00002142991819367
+set pp_ch 0.99910795390933105
+
+restore [locate_data_file bug23706_c2.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp
+if { $pp != $pp_ch } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_22 b/tests/bugs/modalg_5/bug23706_22
new file mode 100644 (file)
index 0000000..b95aa01
--- /dev/null
@@ -0,0 +1,22 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.99999991301930024
+set y 5.00000000157289337
+set z 5.00000535747954842
+set pp_ch 0.99955486819834238
+
+restore [locate_data_file bug23706_c2.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp
+if { $pp != $pp_ch } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_24 b/tests/bugs/modalg_5/bug23706_24
new file mode 100644 (file)
index 0000000..8a2672b
--- /dev/null
@@ -0,0 +1,22 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.99999999837571056
+set y 5.0000000000293724
+set z 5.0000001000463034
+set pp_ch 0.9999392756740122
+
+restore [locate_data_file bug23706_c2.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp
+if { $pp != $pp_ch } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_26 b/tests/bugs/modalg_5/bug23706_26
new file mode 100644 (file)
index 0000000..bed1d5b
--- /dev/null
@@ -0,0 +1,25 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 11.0
+set y -6.0
+set z 5.0
+set pp_ch1 0.22894170490369878
+set pp_ch2 1.7710582950963012
+
+restore [locate_data_file bug23706_c03.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
+
diff --git a/tests/bugs/modalg_5/bug23706_27 b/tests/bugs/modalg_5/bug23706_27
new file mode 100644 (file)
index 0000000..506d935
--- /dev/null
@@ -0,0 +1,24 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.0
+set y 6.0
+set z -3.0
+set pp_ch1 1
+set pp_ch2 1
+
+restore [locate_data_file bug23706_c03.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_28 b/tests/bugs/modalg_5/bug23706_28
new file mode 100644 (file)
index 0000000..05e7454
--- /dev/null
@@ -0,0 +1,26 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 5.0
+set y 8.0
+set z -2.0
+set pp_ch1 1
+set pp_ch2 1
+set pp_ch3 1.1865241781930462
+
+restore [locate_data_file bug23706_c03.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_29 b/tests/bugs/modalg_5/bug23706_29
new file mode 100644 (file)
index 0000000..2d1bc73
--- /dev/null
@@ -0,0 +1,26 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x -4.0
+set y 4.0
+set z 1.0
+set pp_ch1 0
+set pp_ch2 1
+set pp_ch3 1
+
+restore [locate_data_file bug23706_c03.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_3 b/tests/bugs/modalg_5/bug23706_3
new file mode 100755 (executable)
index 0000000..71c5592
--- /dev/null
@@ -0,0 +1,23 @@
+puts "============"
+puts "OCC23706"
+puts "============"
+puts ""
+#########################################################################
+# Cannot project point on curve
+#########################################################################
+
+restore [locate_data_file bug23706_c2.draw] c
+mkedge e1 c
+wire w e1
+polyline pl 3.9 4.9 5 4 5.1 5 4.1 3.9 5 3.9 4.9 5
+mkplane pl pl
+pipe result w pl
+don result
+fit
+checkshape result
+
+set square 1.86489
+
+
+
+
diff --git a/tests/bugs/modalg_5/bug23706_31 b/tests/bugs/modalg_5/bug23706_31
new file mode 100644 (file)
index 0000000..6f53a40
--- /dev/null
@@ -0,0 +1,25 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 11.0
+set y -6.0
+set z 5.0
+set pp_ch1 0.22894170490369881
+set pp_ch2 1.7710582950963012
+
+restore [locate_data_file bug23706_c04.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
+
diff --git a/tests/bugs/modalg_5/bug23706_32 b/tests/bugs/modalg_5/bug23706_32
new file mode 100644 (file)
index 0000000..6f7bda2
--- /dev/null
@@ -0,0 +1,24 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.0
+set y 6.0
+set z -3.0
+set pp_ch1 1
+set pp_ch2 1
+
+restore [locate_data_file bug23706_c04.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_33 b/tests/bugs/modalg_5/bug23706_33
new file mode 100644 (file)
index 0000000..2ca8763
--- /dev/null
@@ -0,0 +1,26 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 5.0
+set y 8.0
+set z -2.0
+set pp_ch1 0.81347582180695399
+set pp_ch2 1
+set pp_ch3 1
+
+restore [locate_data_file bug23706_c04.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_34 b/tests/bugs/modalg_5/bug23706_34
new file mode 100644 (file)
index 0000000..b14a05f
--- /dev/null
@@ -0,0 +1,26 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x -4.0
+set y 4.0
+set z 1.0
+set pp_ch1 1
+set pp_ch2 1
+set pp_ch3 2
+
+restore [locate_data_file bug23706_c04.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_36 b/tests/bugs/modalg_5/bug23706_36
new file mode 100644 (file)
index 0000000..22bc5dd
--- /dev/null
@@ -0,0 +1,25 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 11.0
+set y -6.0
+set z 5.0
+set pp_ch1 0.22894170490369881
+set pp_ch2 1.7205732840814361
+
+restore [locate_data_file bug23706_c05.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
+
diff --git a/tests/bugs/modalg_5/bug23706_37 b/tests/bugs/modalg_5/bug23706_37
new file mode 100644 (file)
index 0000000..4ffd361
--- /dev/null
@@ -0,0 +1,24 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.0
+set y 6.0
+set z -3.0
+set pp_ch1 1
+set pp_ch2 1
+
+restore [locate_data_file bug23706_c05.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_38 b/tests/bugs/modalg_5/bug23706_38
new file mode 100644 (file)
index 0000000..d55ae36
--- /dev/null
@@ -0,0 +1,26 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 5.0
+set y 8.0
+set z -2.0
+set pp_ch1 1
+set pp_ch2 1
+set pp_ch3 1.0371228345434986
+
+restore [locate_data_file bug23706_c05.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_39 b/tests/bugs/modalg_5/bug23706_39
new file mode 100644 (file)
index 0000000..bd679f6
--- /dev/null
@@ -0,0 +1,26 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x -4.0
+set y 4.0
+set z 1.0
+set pp_ch1 0
+set pp_ch2 1
+set pp_ch3 1
+
+restore [locate_data_file bug23706_c05.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_4 b/tests/bugs/modalg_5/bug23706_4
new file mode 100755 (executable)
index 0000000..e84fa2b
--- /dev/null
@@ -0,0 +1,37 @@
+puts "============"
+puts "OCC23706"
+puts "============"
+puts ""
+#########################################################################
+# Cannot project point on curve
+#########################################################################
+
+2dbsplinecurve c1 2 5 0 3 0.2 1 0.3 1 0.4 1 0.5 3 2 0 1 3 -1 1 5 5 1 5 5 1 6 8 1 4 7 1
+2dbsplinecurve c3 2 4 1 3 2 1 3 1 5 3 6 3 1 5.001 5.01 1 5.001 5.01 1 3 9 1 2 11 1
+set info [2dintersect c1 c3]
+
+if { [regexp "Intersection point 1" $info] != 1  } {
+   puts "Error : Intersection should have two points"
+} else {
+   regexp {Intersection point 1 +: +([-0-9.+eE]+)} $info full p11t
+   regexp {Intersection point 1 +: +[-0-9.+eE]+ +([-0-9.+eE]+)} $info full p12t
+}
+
+if { [regexp "Intersection point 2" $info] != 1  } {
+   puts "Error : Intersection should have two points"
+} else {
+   regexp {Intersection point 2 +: +([-0-9.+eE]+)} $info full p21t
+   regexp {Intersection point 2 +: +[-0-9.+eE]+ +([-0-9.+eE]+)} $info full p22t
+}
+
+set p11 [expr round($p11t*10000)]
+set p12 [expr round($p12t*10000)]
+set p21 [expr round($p21t*10000)]
+set p22 [expr round($p22t*10000)]
+
+if { ${p11} != 50024 || ${p12} != 50072 || ${p21} != 40024 || ${p22} != 70012 } {
+    puts "Error : Points of intersection have wrong coordinates"
+} else {
+    puts "OK: Points of intersection are right"
+}
+
diff --git a/tests/bugs/modalg_5/bug23706_40 b/tests/bugs/modalg_5/bug23706_40
new file mode 100644 (file)
index 0000000..00fe14a
--- /dev/null
@@ -0,0 +1,22 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x -3.0
+set y 15.0
+set z -9.0
+set pp_ch 0.31967360381308058
+
+restore [locate_data_file bug23706_c07.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp
+if { $pp != $pp_ch } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_41 b/tests/bugs/modalg_5/bug23706_41
new file mode 100644 (file)
index 0000000..73d21ae
--- /dev/null
@@ -0,0 +1,23 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 11.0
+set y -6.0
+set z 5.0
+set pp_ch 1.7205732840814361
+
+restore [locate_data_file bug23706_c07.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp
+if { $pp != $pp_ch } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
+
diff --git a/tests/bugs/modalg_5/bug23706_42 b/tests/bugs/modalg_5/bug23706_42
new file mode 100644 (file)
index 0000000..7d133ce
--- /dev/null
@@ -0,0 +1,24 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.0
+set y 6.0
+set z -3.0
+set pp_ch1 1
+set pp_ch2 1
+
+restore [locate_data_file bug23706_c07.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_43 b/tests/bugs/modalg_5/bug23706_43
new file mode 100644 (file)
index 0000000..bc482f1
--- /dev/null
@@ -0,0 +1,26 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 5.0
+set y 8.0
+set z -2.0
+set pp_ch1 1
+set pp_ch2 1
+set pp_ch3 1.0371228345434986
+
+restore [locate_data_file bug23706_c07.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_44 b/tests/bugs/modalg_5/bug23706_44
new file mode 100644 (file)
index 0000000..749b807
--- /dev/null
@@ -0,0 +1,26 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x -4.0
+set y 4.0
+set z 1.0
+set pp_ch1 0.034819847916144751
+set pp_ch2 1
+set pp_ch3 1
+
+restore [locate_data_file bug23706_c07.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_45 b/tests/bugs/modalg_5/bug23706_45
new file mode 100644 (file)
index 0000000..9020822
--- /dev/null
@@ -0,0 +1,22 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x -3.0
+set y 15.0
+set z -9.0
+set pp_ch 0.51963826852813999
+
+restore [locate_data_file bug23706_c08.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp
+if { $pp != $pp_ch } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_46 b/tests/bugs/modalg_5/bug23706_46
new file mode 100644 (file)
index 0000000..d5166af
--- /dev/null
@@ -0,0 +1,24 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 11.0
+set y -6.0
+set z 5.0
+set pp_ch1 0.48292970726418566
+set pp_ch2 1.7205732840814361
+
+restore [locate_data_file bug23706_c08.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_47 b/tests/bugs/modalg_5/bug23706_47
new file mode 100644 (file)
index 0000000..abaa114
--- /dev/null
@@ -0,0 +1,24 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.0
+set y 6.0
+set z -3.0
+set pp_ch1 1
+set pp_ch2 1
+
+restore [locate_data_file bug23706_c08.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_48 b/tests/bugs/modalg_5/bug23706_48
new file mode 100644 (file)
index 0000000..8bf2f89
--- /dev/null
@@ -0,0 +1,26 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 5.0
+set y 8.0
+set z -2.0
+set pp_ch1 1
+set pp_ch2 1
+set pp_ch3 1.0371228345434986
+
+restore [locate_data_file bug23706_c08.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_49 b/tests/bugs/modalg_5/bug23706_49
new file mode 100644 (file)
index 0000000..cd75d4d
--- /dev/null
@@ -0,0 +1,26 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x -4.0
+set y 4.0
+set z 1.0
+set pp_ch1 0.087689905182099182
+set pp_ch2 1
+set pp_ch3 1
+
+restore [locate_data_file bug23706_c08.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_5 b/tests/bugs/modalg_5/bug23706_5
new file mode 100755 (executable)
index 0000000..19d169f
--- /dev/null
@@ -0,0 +1,34 @@
+puts "============"
+puts "OCC23706"
+puts "============"
+puts ""
+#########################################################################
+# Cannot project point on curve
+#########################################################################
+
+2dbsplinecurve c1 2 5 0 3 0.2 1 0.3 1 0.4 1 0.5 3 2 0 1 3 -1 1 5 5 1 5 5 1 6 8 1 4 7 1
+2dbsplinecurve c2 2 4 1 3 2 1 3 1 5 3 6 3 1 5 5 1 5 5 1 3 9 1 2 11 1
+set info [2dintersect c1 c2]
+
+if { [regexp "Intersection point 1" $info] != 1  } {
+   puts "Error : Intersection should have two points"
+} else {
+   regexp {Intersection point 1 +: +([-0-9.+eE]+)} $info full p11
+   regexp {Intersection point 1 +: +[-0-9.+eE]+ +([-0-9.+eE]+)} $info full p12
+}
+
+if { [regexp "Intersection point 2" $info] != 1  } {
+   puts "Error : Intersection should have two points"
+} else {
+   regexp {Intersection point 2 +: +([-0-9.+eE]+)} $info full p21t
+   regexp {Intersection point 2 +: +[-0-9.+eE]+ +([-0-9.+eE]+)} $info full p22t
+}
+
+set p21 [expr int($p21t)]
+set p22 [expr int($p22t)]
+
+if { ${p11} != 5 || ${p12} != 5 || ${p21} != 4 || ${p22} != 7 } {
+    puts "Error : Points of intersection have wrong coordinates"
+} else {
+    puts "OK: Points of intersection are right"
+}
diff --git a/tests/bugs/modalg_5/bug23706_50 b/tests/bugs/modalg_5/bug23706_50
new file mode 100644 (file)
index 0000000..18eedbb
--- /dev/null
@@ -0,0 +1,40 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.0
+set y 3.0
+set z 2.0
+set pp_ch1 2.261838779028444
+set pp_ch2 2.7514388736312116
+set pp_ch3 3.5195936992321921
+set pp_ch4 3.9600115496393977
+set pp_ch5 5.4999999987220543
+set pp_ch6 6.8388132447593541
+set pp_ch7 7.8261046366621292
+
+restore [locate_data_file bug23706_c11.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4
+regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5
+regexp {parameter 6 += +([-0-9.+eE]+)} $info full pp6
+regexp {parameter 7 += +([-0-9.+eE]+)} $info full pp7
+if { $pp1 != $pp_ch1 ||
+     $pp2 != $pp_ch2 ||
+     $pp3 != $pp_ch3 ||
+     $pp4 != $pp_ch4 ||
+     $pp5 != $pp_ch5 ||
+     $pp6 != $pp_ch6 ||
+     $pp7 != $pp_ch7} {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_51 b/tests/bugs/modalg_5/bug23706_51
new file mode 100644 (file)
index 0000000..1c2938b
--- /dev/null
@@ -0,0 +1,34 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 5.0
+set y 7.0
+set z 8.0
+set pp_ch1 2.8126840147763663
+set pp_ch2 3.5195936992321926
+set pp_ch3 3.9600115496393977
+set pp_ch4 5.4999999987220543
+set pp_ch5 7.2883607799598096
+set pp_ch6 1
+
+restore [locate_data_file bug23706_c11.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4
+regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5
+regexp {parameter 6 += +([-0-9.+eE]+)} $info full pp6
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 ||
+     $pp3 != $pp_ch3 || $pp4 != $pp_ch4 ||
+     $pp5 != $pp_ch5 || $pp6 != $pp_ch6} {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_52 b/tests/bugs/modalg_5/bug23706_52
new file mode 100644 (file)
index 0000000..fd19bca
--- /dev/null
@@ -0,0 +1,30 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 11.0
+set y -2.0
+set z -2.0
+set pp_ch1 2.9473269594602054
+set pp_ch2 4.4416670680933228
+set pp_ch3 5.4999999987220543
+set pp_ch4 6.6582576262308306
+set pp_ch5 7.7414573419084736
+
+restore [locate_data_file bug23706_c11.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4
+regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 || $pp4 != $pp_ch4 || $pp5 != $pp_ch5} {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_53 b/tests/bugs/modalg_5/bug23706_53
new file mode 100644 (file)
index 0000000..4fe4224
--- /dev/null
@@ -0,0 +1,40 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.0
+set y 3.0
+set z 2.0
+set pp_ch1 1.1738953633378706
+set pp_ch2 2.1611867552406454
+set pp_ch3 3.5000000012779413
+set pp_ch4 5.0399884503606023
+set pp_ch5 5.4804063007678074
+set pp_ch6 6.2485611263687888
+set pp_ch7 6.7381612209715556
+
+restore [locate_data_file bug23706_c12.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4
+regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5
+regexp {parameter 6 += +([-0-9.+eE]+)} $info full pp6
+regexp {parameter 7 += +([-0-9.+eE]+)} $info full pp7
+if { $pp1 != $pp_ch1 ||
+     $pp2 != $pp_ch2 ||
+     $pp3 != $pp_ch3 ||
+     $pp4 != $pp_ch4 ||
+     $pp5 != $pp_ch5 ||
+     $pp6 != $pp_ch6 ||
+     $pp7 != $pp_ch7} {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_54 b/tests/bugs/modalg_5/bug23706_54
new file mode 100644 (file)
index 0000000..e3cbb58
--- /dev/null
@@ -0,0 +1,34 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 5.0
+set y 7.0
+set z 8.0
+set pp_ch1 1.7116392200401909
+set pp_ch2 3.5000000012779413
+set pp_ch3 5.0399884503606023
+set pp_ch4 5.4804063007678074
+set pp_ch5 6.1873159852236332
+set pp_ch6 8
+
+restore [locate_data_file bug23706_c12.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4
+regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5
+regexp {parameter 6 += +([-0-9.+eE]+)} $info full pp6
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 ||
+     $pp3 != $pp_ch3 || $pp4 != $pp_ch4 ||
+     $pp5 != $pp_ch5 || $pp6 != $pp_ch6} {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_55 b/tests/bugs/modalg_5/bug23706_55
new file mode 100644 (file)
index 0000000..833dd50
--- /dev/null
@@ -0,0 +1,30 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 11.0
+set y -2.0
+set z -2.0
+set pp_ch1 1.2585426580915264
+set pp_ch2 2.3417423737691694
+set pp_ch3 3.499999996505935
+set pp_ch4 4.5583329319066772
+set pp_ch5 6.052673040539795
+
+restore [locate_data_file bug23706_c12.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4
+regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 || $pp4 != $pp_ch4 || $pp5 != $pp_ch5} {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_56 b/tests/bugs/modalg_5/bug23706_56
new file mode 100644 (file)
index 0000000..585a783
--- /dev/null
@@ -0,0 +1,34 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.0
+set y 3.0
+set z 2.0
+set pp_ch1 1.8318851868378956
+set pp_ch2 3.0397214383562297
+set pp_ch3 5.5
+set pp_ch4 6.8388132447593541
+set pp_ch5 7.8261046366621292
+
+restore [locate_data_file bug23706_c13.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4
+regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5
+if { $pp1 != $pp_ch1 ||
+     $pp2 != $pp_ch2 ||
+     $pp3 != $pp_ch3 ||
+     $pp4 != $pp_ch4 ||
+     $pp5 != $pp_ch5} {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_57 b/tests/bugs/modalg_5/bug23706_57
new file mode 100644 (file)
index 0000000..d075ef7
--- /dev/null
@@ -0,0 +1,29 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 5.0
+set y 7.0
+set z 8.0
+set pp_ch1 3.0397214383562297
+set pp_ch2 5.5
+set pp_ch3 7.2883607799598096
+set pp_ch4 1
+
+restore [locate_data_file bug23706_c13.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 ||
+     $pp3 != $pp_ch3 || $pp4 != $pp_ch4} {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_58 b/tests/bugs/modalg_5/bug23706_58
new file mode 100644 (file)
index 0000000..2be2338
--- /dev/null
@@ -0,0 +1,30 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 11.0
+set y -2.0
+set z -2.0
+set pp_ch1 2.2389225099869194
+set pp_ch2 3.219764556283669
+set pp_ch3 5.5
+set pp_ch4 6.6582576262308306
+set pp_ch5 7.7414573419084736
+
+restore [locate_data_file bug23706_c13.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4
+regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 || $pp4 != $pp_ch4 || $pp5 != $pp_ch5} {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_59 b/tests/bugs/modalg_5/bug23706_59
new file mode 100644 (file)
index 0000000..58ded81
--- /dev/null
@@ -0,0 +1,34 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 3.0
+set y 3.0
+set z 2.0
+set pp_ch1 1.1738953633378706
+set pp_ch2 2.1611867552406454
+set pp_ch3 3.5000000000000004
+set pp_ch4 5.9602785616437703
+set pp_ch5 7.1681148131621049
+
+restore [locate_data_file bug23706_c14.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4
+regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5
+if { $pp1 != $pp_ch1 ||
+     $pp2 != $pp_ch2 ||
+     $pp3 != $pp_ch3 ||
+     $pp4 != $pp_ch4 ||
+     $pp5 != $pp_ch5} {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_6 b/tests/bugs/modalg_5/bug23706_6
new file mode 100755 (executable)
index 0000000..8314b56
--- /dev/null
@@ -0,0 +1,18 @@
+puts "============"
+puts "OCC23706"
+puts "============"
+puts ""
+#########################################################################
+# Cannot project point on curve
+#########################################################################
+
+pload XSDRAW
+2dbsplinecurve cc 3 2 0 4 1 4 -1 -1 1 0 -1 1 0 0 1 0 0 1
+offset2dcurve o cc .5
+set info [length o]
+regexp {The length o is+ +([-0-9.+eE]+)} $info full ll
+set ll_check 2.3717833300483151
+
+if { [expr 1.*abs($ll_check - $ll)/$ll_check] > 0.01 } {
+   puts "Error : The lenght of result shape is $ll"
+}
diff --git a/tests/bugs/modalg_5/bug23706_60 b/tests/bugs/modalg_5/bug23706_60
new file mode 100644 (file)
index 0000000..5ec0c89
--- /dev/null
@@ -0,0 +1,29 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 5.0
+set y 7.0
+set z 8.0
+set pp_ch1 1.7116392200401909
+set pp_ch2 3.5000000000000004
+set pp_ch3 5.9602785616437703
+set pp_ch4 8
+
+restore [locate_data_file bug23706_c14.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 ||
+     $pp3 != $pp_ch3 || $pp4 != $pp_ch4} {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_61 b/tests/bugs/modalg_5/bug23706_61
new file mode 100644 (file)
index 0000000..8057244
--- /dev/null
@@ -0,0 +1,30 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+####################################
+## Cannot project point on curve
+####################################
+
+set x 11.0
+set y -2.0
+set z -2.0
+set pp_ch1 1.2585426580915264
+set pp_ch2 2.3417423737691694
+set pp_ch3 3.5000000000000004
+set pp_ch4 5.7802354437163306
+set pp_ch5 6.761077490013081
+
+restore [locate_data_file bug23706_c14.draw] c
+set info [proj c $x $y $z]
+
+regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1
+regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2
+regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3
+regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4
+regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5
+if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 || $pp4 != $pp_ch4 || $pp5 != $pp_ch5} {
+   puts "Error : Projection is not correct"
+} else {
+   puts "OK: Projection is correct"
+}
diff --git a/tests/bugs/modalg_5/bug23706_7 b/tests/bugs/modalg_5/bug23706_7
new file mode 100755 (executable)
index 0000000..ae4941f
--- /dev/null
@@ -0,0 +1,19 @@
+puts "============"
+puts "OCC23706"
+puts "============"
+puts ""
+#########################################################################
+# Cannot project point on curve
+#########################################################################
+
+pload XSDRAW
+bsplinecurve cc 3 2 0 4 1 4 -1 -1 2 1 0 -1 2 1 0 0 2 1 0 0 2 1
+point pp 0 0 1
+offsetcurve o cc .5 pp
+set info [length o]
+regexp {The length o is+ +([-0-9.+eE]+)} $info full ll
+set ll_check 2.3717833300483151
+
+if { [expr 1.*abs($ll_check - $ll)/$ll_check] > 0.01 } {
+   puts "Error : The lenght of result shape is $ll"
+}
diff --git a/tests/bugs/modalg_5/bug23706_8 b/tests/bugs/modalg_5/bug23706_8
new file mode 100755 (executable)
index 0000000..2f7a3c6
--- /dev/null
@@ -0,0 +1,25 @@
+puts "============"
+puts "OCC23706"
+puts "============"
+puts ""
+#########################################################################
+# Cannot project point on curve
+#########################################################################
+
+bsplinecurve r2 2 4 1 3 2 2 3 1 4 3 0 8 0 1 2 8 2 1 4 8 3 1 4 8 3 1 6 8 4 1 10 8 10 1
+mkedge spine r2
+wire spine spine
+circle profile 0 8 0 1 0 1 0.51
+mkedge profile profile
+wire profile profile
+mkplane profile profile
+pipe p spine profile
+explode p f
+mksurface ss p_2
+bsplinecurve r3 2 5 1 3 2 1 3 1 4 1 5 3 9 7 3 1 5.447213595499958 9 2.105572809000084 1 4.223606797749979 8 2.552786404500042 1 4.223606797749979 8 2.552786404500042 1 3 7 3 1 2 5 3 1
+set info [intersect res r3 ss]
+if { [regexp "res_1" $info] != 1 || [regexp "res_2" $info] != 1 || [regexp "res_3" $info] != 0 } {
+    puts "Error : Wrong number of points of intersection. Should contain two points"
+} else {
+    puts "OK: Number of points of intersection is valid"
+}
diff --git a/tests/bugs/modalg_5/bug23706_9 b/tests/bugs/modalg_5/bug23706_9
new file mode 100755 (executable)
index 0000000..bb7d45d
--- /dev/null
@@ -0,0 +1,19 @@
+puts "============"
+puts "OCC23706"
+puts "============"
+puts ""
+#########################################################################
+# Cannot project point on curve
+#########################################################################
+
+pload XSDRAW
+2dbsplinecurve c1 2 5 0 3 0.2 1 0.3 1 0.4 1 0.5 3 2 0 1 3 -1 1 5 5 1 5 5 1 6 8 1 4 7 1
+offset2dcurve o1 c1 2
+
+set info [length o1]
+regexp {The length o1 is+ +([-0-9.+eE]+)} $info full ll
+set ll_check 19.244437838214424
+
+if { [expr 1.*abs($ll_check - $ll)/$ll_check] > 0.01 } {
+   puts "Error : The lenght of result shape is $ll"
+}
diff --git a/tests/bugs/moddata_3/bug23706 b/tests/bugs/moddata_3/bug23706
new file mode 100644 (file)
index 0000000..4209273
--- /dev/null
@@ -0,0 +1,20 @@
+puts "========"
+puts "OCC23706"
+puts "========"
+puts ""
+###############################
+## Cannot project point 3.9999965207720098 5.0000000629157348 5.0002142991819367 on curve from attached file.
+###############################
+
+set BugNumber OCC23706
+
+restore [locate_data_file bug23706_c.draw] c 
+proj c 3.9999965207720098 5.0000000629157348 5.0002142991819367
+
+if {[isdraw ext_1] == 1} {
+  dump ext_1
+  puts "${BugNumber} OK"
+} else {
+  puts "Faulty ${BugNumber}"
+}
\ No newline at end of file