0024694: Wrong processing of two surfaces (implicit and parametric) having tangential intersection: it is not found by intersector

src/Contap/Contap_Contour.cxx

@@ -1491,7 +1491,8 @@ void Contap_Contour::Perform
   if (seqpdep.Length() != 0 || seqpins.Length() != 0) {
 
     Contap_TheIWalking iwalk(Preci,Fleche,Pas);
-    iwalk.Perform(seqpdep,seqpins,mySFunc ,Surf);
+    IntSurf_SequenceOfInteriorPoint seqptang; //dummy
+    iwalk.Perform(seqpdep,seqpins,seqptang,mySFunc ,Surf);
     if(!iwalk.IsDone()) {
       return;
     }

src/Contap/Contap_SurfFunction.cdl

@@ -148,6 +148,12 @@ is
     	returns Boolean from Standard 
     	is static;
     
+    IsTangentSmooth(me : in out)
+
+    	returns Boolean from Standard 
+	---C++: inline
+    	is static;
+    
 
     Direction3d(me: in out)
 
@@ -167,6 +173,32 @@ is
     	is static;
 
 
+    DerivativesAndNormalOnPSurf(me: in out; D1U :    out Vec from gp;
+    	    	    	    	    D1V :    out Vec from gp;
+				    Normal : out Vec from gp;
+    	    	    	    	    D2U :    out Vec from gp;
+    	    	    	    	    D2V :    out Vec from gp;
+    	    	    	    	    D2UV :   out Vec from gp)
+    	returns Boolean from Standard
+	---C++: inline
+    	is static;
+    
+    DerivativesAndNormalOnISurf(me; D1U :    out Vec from gp;
+    	    	    	    	    D1V :    out Vec from gp;
+				    Normal : out Vec from gp;
+    	    	    	    	    D2U :    out Vec from gp;
+    	    	    	    	    D2V :    out Vec from gp;
+    	    	    	    	    D2UV :   out Vec from gp)
+    	returns Boolean from Standard
+	---C++: inline
+    	is static;
+    
+    SquareTangentError(me)
+    	returns Real from Standard
+	---C++: inline
+	is static;
+
+
     FunctionType(me)
     
     	returns TFunction from Contap

src/Contap/Contap_SurfFunction.lxx

@@ -100,6 +100,36 @@ inline Standard_Real Contap_SurfFunction::Angle () const
   return myAng;
 }
 
+inline Standard_Boolean Contap_SurfFunction::IsTangentSmooth()
+{
+  return Standard_False;
+}
+
+inline Standard_Boolean Contap_SurfFunction::DerivativesAndNormalOnPSurf(gp_Vec&,
+                                                                          gp_Vec&,
+                                                                          gp_Vec&,
+                                                                          gp_Vec&,
+                                                                          gp_Vec&,
+                                                                          gp_Vec&)
+{
+  return Standard_False;
+}
+
+inline Standard_Boolean Contap_SurfFunction::DerivativesAndNormalOnISurf(gp_Vec&,
+                                                                         gp_Vec&,
+                                                                         gp_Vec&,
+                                                                         gp_Vec&,
+                                                                         gp_Vec&,
+                                                                         gp_Vec&) const
+{
+  return Standard_False;
+}
+
+inline Standard_Real Contap_SurfFunction::SquareTangentError() const
+{
+  return 0.;
+}
+
 inline Contap_TFunction Contap_SurfFunction::FunctionType () const
 {
   return myType;

src/IntImp/IntImp_ZerImpFunc.cdl

@@ -141,6 +141,12 @@ is
 
     	is static;
     
+    IsTangentSmooth(me : in out)
+
+    	returns Boolean from Standard 
+
+    	is static;
+    
 
     Direction3d(me: in out)
 
@@ -161,6 +167,29 @@ is
     	raises UndefinedDerivative from StdFail
     	is static;
     
+    DerivativesAndNormalOnPSurf(me: in out; D1U :    out Vec from gp;
+    	    	    	    	    D1V :    out Vec from gp;
+				    Normal : out Vec from gp;
+    	    	    	    	    D2U :    out Vec from gp;
+    	    	    	    	    D2V :    out Vec from gp;
+    	    	    	    	    D2UV :   out Vec from gp)
+    	returns Boolean from Standard
+    	is static;
+    
+    DerivativesAndNormalOnISurf(me; D1U :    out Vec from gp;
+    	    	    	    	    D1V :    out Vec from gp;
+				    Normal : out Vec from gp;
+    	    	    	    	    D2U :    out Vec from gp;
+    	    	    	    	    D2V :    out Vec from gp;
+    	    	    	    	    D2UV :   out Vec from gp)
+    	returns Boolean from Standard
+    	is static;
+    
+    SquareTangentError(me)
+    	returns Real from Standard
+	---C++: inline
+
+	is static;
 
     PSurface(me)
 
@@ -170,7 +199,6 @@ is
 
     	is static;
     
-
     ISurface(me)
 
     	returns TheISurface
@@ -195,6 +223,8 @@ fields
   tgdu        : Real    from Standard;
   tgdv        : Real    from Standard;
   gradient    : Vec     from gp;
+  d1u_isurf   : Vec     from gp;
+  d1v_isurf   : Vec     from gp;
   derived     : Boolean from Standard;
   d1u         : Vec     from gp;
   d1v         : Vec     from gp;

src/IntImp/IntImp_ZerImpFunc.gxx

@@ -92,7 +92,7 @@ Standard_Boolean IntImp_ZerImpFunc::Values(const math_Vector& X,
   v = X(2);
   ThePSurfaceTool::D1(SURF, u, v, pntsol, d1u, d1v);
   TheISurfaceTool::ValueAndGradient(FUNC, pntsol.X(), pntsol.Y(), pntsol.Z(), 
-                                    valf, gradient);
+                                    valf, gradient, d1u_isurf, d1v_isurf);
   F(1) = valf;
   D(1,1) = d1u.Dot(gradient);
   D(1,2) = d1v.Dot(gradient);
@@ -118,17 +118,133 @@ Standard_Boolean IntImp_ZerImpFunc::IsTangent()
     Standard_Real N2d1v  = d1v.SquareMagnitude();
     tangent =(tgdu * tgdu <= N2grad_EpsAng2 * N2d1v) && 
              (tgdv * tgdv <= N2grad_EpsAng2 * N2d1u);
+
+#ifdef DRAW
+    if (tangent)
+    {
+      Standard_Real SumSquare = tgdu * tgdu + tgdv * tgdv;
+      gp_Vec NormSURF = d1u ^ d1v;
+      cout<<"Tangent !!!"<<endl;
+      cout<<"SumSquare = "<<SumSquare<<endl;
+      cout<<"u = "<<u<<", v = "<<v<<endl;
+      cout<<"pnt = "<<pntsol.X()<<" "<<pntsol.Y()<<" "<<pntsol.Z()<<endl;
+    }
+#endif
+    
     if(!tangent) {
       d3d.SetLinearForm(tgdu,d1u,tgdv,d1v);
       d2d = gp_Dir2d(tgdu, tgdv);
       if (d3d.Magnitude() <= Tolpetit) {    // jag
-        tangent = Standard_True;
+	tangent = Standard_True;
       }
-    }    
+      
+#ifdef DRAW
+      gp_Vec NormSURF = d1u ^ d1v;
+      cout<<"u = "<<u<<", v = "<<v<<endl;
+      cout<<"pnt = "<<pntsol.X()<<" "<<pntsol.Y()<<" "<<pntsol.Z()<<endl;
+      cout<<"d3d = "<<d3d.X()<<" "<<d3d.Y()<<" "<<d3d.Z()<<endl;
+#endif
+    }
   }
   return tangent;
 }
 
+Standard_Boolean IntImp_ZerImpFunc::IsTangentSmooth()
+{
+  const Standard_Real TolCrossProd = 1.e-15;
+  
+  if (!computed) {
+    computed = Standard_True;
+    if(!derived) {
+      ThePSurfaceTool::D1(SURF, u, v, pntsol, d1u, d1v);
+      derived = Standard_True;
+    }
+
+    tgdu =  gradient.Dot(d1v);
+    tgdv = -gradient.Dot(d1u);
+    
+    gp_Vec NormPSurf = d1u ^ d1v;
+    gp_Vec NormISurf = d1u_isurf ^ d1v_isurf;
+    NormPSurf.Normalize();
+    NormISurf.Normalize();
+#ifdef DRAW
+    Standard_Real theAngle = NormPSurf.Angle(NormISurf);
+#endif
+    gp_Vec CrossProd = NormPSurf ^ NormISurf;
+    Standard_Real CrossProdSqMagn = CrossProd.SquareMagnitude();
+    
+#ifdef DRAW
+    cout<<"Angle = "<<theAngle<<endl;
+    cout<<"CrossProdSqMagn = "<<CrossProdSqMagn<<endl;
+#endif
+
+    tangent = (CrossProdSqMagn <= TolCrossProd);
+    
+#ifdef DRAW
+    if (tangent)
+    {
+      Standard_Real SumSquare = tgdu * tgdu + tgdv * tgdv;
+      gp_Vec NormSURF = d1u ^ d1v;
+      cout<<"Tangent !!!"<<endl;
+      cout<<"SumSquare = "<<SumSquare<<endl;
+      cout<<"u = "<<u<<", v = "<<v<<endl;
+      cout<<"pnt = "<<pntsol.X()<<" "<<pntsol.Y()<<" "<<pntsol.Z()<<endl;
+    }
+#endif
+    
+    if(!tangent) {
+      d3d.SetLinearForm(tgdu,d1u,tgdv,d1v);
+      d2d = gp_Dir2d(tgdu, tgdv);
+      if (d3d.Magnitude() <= Tolpetit) {    // jag
+	tangent = Standard_True;
+      }
+#ifdef DRAW
+      gp_Vec NormSURF = d1u ^ d1v;
+      cout<<"u = "<<u<<", v = "<<v<<endl;
+      cout<<"pnt = "<<pntsol.X()<<" "<<pntsol.Y()<<" "<<pntsol.Z()<<endl;
+      cout<<"d3d = "<<d3d.X()<<" "<<d3d.Y()<<" "<<d3d.Z()<<endl;
+#endif
+    }
+  }
+  return tangent;
+}
+
+Standard_Boolean IntImp_ZerImpFunc::DerivativesAndNormalOnPSurf(gp_Vec& D1U,
+                                                    gp_Vec& D1V,
+                                                    gp_Vec& Normal,
+                                                    gp_Vec& D2U,
+                                                    gp_Vec& D2V,
+                                                    gp_Vec& D2UV)
+{
+  ThePSurfaceTool::D2(SURF, u, v, pntsol, D1U, D1V, D2U, D2V, D2UV);
+  Normal = D1U ^ D1V;
+  Standard_Real theMagnitude = Normal.Magnitude();
+  if (theMagnitude <= gp::Resolution())
+    return Standard_False;
+  else
+    Normal /= theMagnitude;
+
+  return Standard_True;
+}
+
+Standard_Boolean IntImp_ZerImpFunc::DerivativesAndNormalOnISurf(gp_Vec& D1U,
+                                                    gp_Vec& D1V,
+                                                    gp_Vec& Normal,
+                                                    gp_Vec& D2U,
+                                                    gp_Vec& D2V,
+                                                    gp_Vec& D2UV) const
+{
+  TheISurfaceTool::D2(FUNC, pntsol, D1U, D1V, D2U, D2V, D2UV);
+  Normal = D1U ^ D1V;
+  Standard_Real theMagnitude = Normal.Magnitude();
+  if (theMagnitude <= gp::Resolution())
+    return Standard_False;
+  else
+    Normal /= theMagnitude;
+
+  return Standard_True;
+}
+
 #undef EpsAng
 #undef EpsAng2
 #undef Tolpetit

src/IntImp/IntImp_ZerImpFunc.lxx

@@ -56,6 +56,11 @@ inline const gp_Dir2d& IntImp_ZerImpFunc::Direction2d()
   return d2d;
 }
 
+inline Standard_Real IntImp_ZerImpFunc::SquareTangentError() const
+{
+  return (tgdu * tgdu + tgdv * tgdv);
+}
+
 inline const ThePSurface& IntImp_ZerImpFunc::PSurface() const 
 {
   return  (*((ThePSurface *)(surf)));

src/IntPatch/IntPatch_ImpPrmIntersection.cxx

@@ -413,7 +413,7 @@ void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur
 {
   Standard_Boolean reversed, procf, procl, dofirst, dolast;
   Standard_Integer indfirst = 0, indlast = 0, ind2, NbSegm;
-  Standard_Integer NbPointIns, NbPointRst, Nblines, Nbpts, NbPointDep;
+  Standard_Integer NbPointIns, NbPointsTang, NbPointRst, Nblines, Nbpts, NbPointDep;
   Standard_Real U1,V1,U2,V2,paramf,paraml,currentparam;
 
   IntPatch_TheSegmentOfTheSOnBounds thesegm;
@@ -528,6 +528,7 @@ void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur
   //
   IntSurf_SequenceOfPathPoint seqpdep;
   IntSurf_SequenceOfInteriorPoint seqpins;
+  IntSurf_SequenceOfInteriorPoint seqptang;
   //
   NbPointRst = solrst.NbPoints();
   TColStd_Array1OfInteger Destination(1,NbPointRst+1); Destination.Init(0);
@@ -599,6 +600,9 @@ void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur
     for (Standard_Integer i=1; i <= NbPointIns; i++) {
       seqpins.Append(solins.Value(i));
     }
+    NbPointsTang = solins.NbTangentPoints();
+    for (Standard_Integer i = 1; i <= NbPointsTang; i++)
+      seqptang.Append(solins.TangentPoint(i));
   }
   //
   NbPointDep=seqpdep.Length();
@@ -606,10 +610,10 @@ void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_HSurface)& Sur
   if (NbPointDep || NbPointIns) {
     IntPatch_TheIWalking iwalk(TolTang,Fleche,Pas);
     if (!reversed) {
-      iwalk.Perform(seqpdep,seqpins,Func,Surf2);
+      iwalk.Perform(seqpdep,seqpins,seqptang,Func,Surf2);
     }
     else {
-      iwalk.Perform(seqpdep,seqpins,Func,Surf1,Standard_True);
+      iwalk.Perform(seqpdep,seqpins,seqptang,Func,Surf1,Standard_True);
     }
     if(!iwalk.IsDone()) {
       return;

src/IntStart/IntStart_SearchInside.cdl

@@ -68,7 +68,20 @@ is
     NbPoints(me)
     
 	---Purpose: Returns the number of points.
-	--          The exception NotDone if raised if IsDone 
+	--          The exception NotDone is raised if IsDone 
+	--          returns False.
+    
+    	returns Integer
+	---C++: inline
+	
+	raises NotDone from StdFail
+	
+	is static;
+
+    NbTangentPoints(me)
+    
+	---Purpose: Returns the number of tangent points.
+	--          The exception NotDone is raised if IsDone 
 	--          returns False.
     
     	returns Integer
@@ -96,10 +109,29 @@ is
 
     	is static;
 
+    TangentPoint(me; Index: Integer)
+    
+	---Purpose: Returns the tangent point of range Index.
+	--          
+	--          The exception NotDone is raised if IsDone 
+	--          returns False.
+	--          The exception OutOfRange if raised if
+	--          Index <= 0 or Index > NbTangentPoints.
+
+    	returns InteriorPoint from IntSurf
+	---C++: return const&
+	---C++: inline
+	
+	raises NotDone from StdFail,
+	       OutOfRange from Standard
+
+    	is static;
+
 
 fields
 
-    done : Boolean                 from Standard;
+    done            : Boolean                 from Standard;
-    list : SequenceOfInteriorPoint from IntSurf;
+    list            : SequenceOfInteriorPoint from IntSurf;
+    myTangentPoints : SequenceOfInteriorPoint from IntSurf;
 
 end SearchInside;

src/IntStart/IntStart_SearchInside.gxx

@@ -25,6 +25,11 @@
 
 #include <TopAbs_State.hxx>
 
+#ifdef DRAW
+#include <DrawTrSurf.hxx>
+#endif
+
+
 IntStart_SearchInside::IntStart_SearchInside (): done(Standard_False)
 {}
 
@@ -49,6 +54,7 @@ void IntStart_SearchInside::Perform (TheFunction& Func,
 
   done = Standard_False;
   list.Clear();
+  myTangentPoints.Clear();
   Standard_Real aBinf[2], aBsup[2], aUVap[2], atoler[2];
   math_Vector Binf(aBinf,1,2), Bsup(aBsup,1,2), UVap(aUVap,1,2), toler(atoler,1,2);
   gp_Pnt psol;
@@ -118,6 +124,13 @@ void IntStart_SearchInside::Perform (TheFunction& Func,
   //--
 
   
+#ifdef DRAW
+  char* name = new char[100];
+#endif
+
+  gp_Vec NullDir(0.,0.,0.);
+  gp_Vec2d NullDir2d(0.,0.);
+  
   for (i=1; i <= Nbsample+12; i++) {
     gp_Pnt2d s2d;
     gp_Pnt s3d;
@@ -176,9 +189,9 @@ void IntStart_SearchInside::Perform (TheFunction& Func,
       Rsnld.Perform(Func,UVap,Binf,Bsup);
       if (Rsnld.IsDone()) {
 	if (Abs(Func.Root()) <= Tol) {
+          psol = Func.Point();
+          Rsnld.Root(UVap);
 	  if (!Func.IsTangent()) {
-	    psol = Func.Point();
-	    Rsnld.Root(UVap);
 	    // On regarde si le point trouve est bien un nouveau point.	  
 	    j = 1;
 	    nbpt = list.Length();
@@ -207,9 +220,54 @@ void IntStart_SearchInside::Perform (TheFunction& Func,
 		list.Append(IntSurf_InteriorPoint(psol,UVap(1),UVap(2),
 						  Func.Direction3d(),
 						  Func.Direction2d()));
+#ifdef DRAW
+                gp_Pnt2d p2dsol(UVap(1),UVap(2));
+                sprintf(name, "pp%d", list.Length());
+                DrawTrSurf::Set(name, p2dsol);
+                sprintf(name, "p%d", list.Length());
+                DrawTrSurf::Set(name, psol);
+#endif
 	      }
 	    }
-	  }
+	  }//if (!Func.IsTangent())
+          else //tangent point
+          {
+            Standard_Boolean IsNewTangentPoint = Standard_True;
+            for (j = 1; j <= myTangentPoints.Length(); j++)
+            {
+	      const IntSurf_InteriorPoint&  IPj = myTangentPoints(j);
+	      const gp_Pnt& Pj = IPj.Value();
+	      if (   (Abs(Pj.X()-psol.X()) <= Epsilon)
+		  && (Abs(Pj.Y()-psol.Y()) <= Epsilon)
+		  && (Abs(Pj.Z()-psol.Z()) <= Epsilon)
+		  && (Abs(UVap(1)-IPj.UParameter()) <= toler1)
+		  && (Abs(UVap(2)-IPj.VParameter()) <= toler2) )
+              {
+		IsNewTangentPoint = Standard_False;
+                break;
+	      }
+            }
+            
+            if (IsNewTangentPoint)
+            {
+	      situ = T->Classify(gp_Pnt2d(UVap(1),UVap(2)),
+				 Maxtoler1toler2,Standard_False);   //-- ,Standard_False pour ne pas recadrer on Periodic
+	      if (situ == TopAbs_IN) {
+                myTangentPoints.Append(IntSurf_InteriorPoint(psol,UVap(1),UVap(2),
+                                                             NullDir, NullDir2d));
+#ifdef DRAW
+                gp_Pnt2d p2dsol(UVap(1),UVap(2));
+                sprintf(name, "tpp%d", myTangentPoints.Length());
+                //DrawTrSurf::Set(name, p2dsol);
+                sprintf(name, "tp%d", myTangentPoints.Length());
+                //DrawTrSurf::Set(name, psol);
+                cout<<myTangentPoints.Length()<<" tangent point :"<<endl;
+                cout<<"U = "<<UVap(1)<<", V = "<<UVap(2)<<endl;
+                cout<<"Dist = "<<aRootVal<<endl<<endl;
+#endif
+              }
+            }
+          }
 	}
       }
     }
@@ -231,6 +289,7 @@ void IntStart_SearchInside::Perform (TheFunction& Func,
 {
   done = Standard_False;
   list.Clear();
+  myTangentPoints.Clear();
   math_Vector Binf(1,2), Bsup(1,2), toler(1,2);
 
   Binf(1) = ThePSurfaceTool::FirstUParameter(PS);

src/IntStart/IntStart_SearchInside.lxx

@@ -27,6 +27,12 @@ inline Standard_Integer IntStart_SearchInside::NbPoints () const {
   return list.Length();
 }
 
+inline Standard_Integer IntStart_SearchInside::NbTangentPoints () const {
+
+  if (!done) {StdFail_NotDone::Raise();}
+  return myTangentPoints.Length();
+}
+
 
 inline const IntSurf_InteriorPoint& IntStart_SearchInside::Value 
      (const Standard_Integer Index) const
@@ -36,3 +42,10 @@ inline const IntSurf_InteriorPoint& IntStart_SearchInside::Value
   return list.Value(Index);
 }
 
+inline const IntSurf_InteriorPoint& IntStart_SearchInside::TangentPoint
+     (const Standard_Integer Index) const
+
+{
+  if (!done) {StdFail_NotDone::Raise();}
+  return myTangentPoints.Value(Index);
+}

src/IntSurf/IntSurf_PntOn2S.cdl

@@ -22,7 +22,9 @@ class PntOn2S from IntSurf
 	--          parametric coordinates.
 
 
-uses Pnt from gp
+uses
+    Pnt   from gp,
+    Pnt2d from gp
 
 is
 
@@ -91,6 +93,13 @@ is
 
 	is static;
 
+    ValueOnSurface(me; OnFirst: Boolean from Standard)
+    
+	---Purpose: Returns the point in 2d space of one of the surfaces.
+
+    	returns Pnt2d from gp
+	is static;
+
 
 
     ParametersOnS1(me; U1,V1: out Real from Standard)
@@ -110,6 +119,13 @@ is
 
     	is static;
 
+    ParametersOnSurface(me; OnFirst: Boolean from Standard;
+    	    	    	    U,V: out Real from Standard)
+    
+	---Purpose: Returns the parameters of the point in the
+	--          parametric space of one of the surface.
+
+    	is static;
 
     Parameters(me; U1,V1,U2,V2: out Real from Standard)
 
@@ -119,7 +135,7 @@ is
 
     	is static;
 
-    IsSame(me; theOterPoint : PntOn2S from IntSurf;
+    IsSame(me; theOtherPoint : PntOn2S from IntSurf;
             theTol3D, theTol2D: Real from Standard = 0.0)
             
   ---Purpose: Returns TRUE if 2D- and 3D-coordinates of theOterPoint are equal to

src/IntSurf/IntSurf_PntOn2S.cxx

@@ -48,6 +48,29 @@ void IntSurf_PntOn2S::SetValue (const Standard_Boolean OnFirst,
   }
 }
 
+gp_Pnt2d IntSurf_PntOn2S::ValueOnSurface(const Standard_Boolean OnFirst) const
+{
+  gp_Pnt2d PointOnSurf;
+  if (OnFirst)
+    PointOnSurf.SetCoord(u1,v1);
+  else
+    PointOnSurf.SetCoord(u2,v2);
+  return PointOnSurf;
+}
+
+void IntSurf_PntOn2S::ParametersOnSurface(const Standard_Boolean OnFirst,
+                                          Standard_Real& U,
+                                          Standard_Real& V) const
+{
+  if (OnFirst) {
+    U = u1;
+    V = v1;
+  }
+  else {
+    U = u2;
+    V = v2;
+  }
+}
 
 Standard_Boolean IntSurf_PntOn2S::IsSame( const IntSurf_PntOn2S& theOterPoint,
                                           const Standard_Real theTol3D,

src/IntSurf/IntSurf_QuadricTool.cdl

@@ -22,7 +22,8 @@ class QuadricTool from IntSurf
 	--          as implicit surface.
 
 uses Quadric from IntSurf,
-     Vec     from gp
+     Vec     from gp,
+     Pnt     from gp
 
 is
 
@@ -47,11 +48,22 @@ is
 
     ValueAndGradient(myclass; Quad: Quadric from IntSurf;
                      X, Y, Z: Real from Standard;
-                     Val: out Real from Standard; Grad: out Vec from gp);
+                     Val: out Real from Standard;
+    	    	     Grad: out Vec from gp;
+		     D1U_ISurf, D1V_ISurf: out Vec from gp);
 		     
     	---Purpose: Returns the value and the gradient.
 
-	---C++: inline
+    D2(myclass; Quad:  Quadric from IntSurf;
+                Point: Pnt from gp;
+    	    	D1U:   out Vec from gp;
+    	    	D1V:   out Vec from gp;
+    	    	D2U:   out Vec from gp;
+    	    	D2V:   out Vec from gp;
+    	    	D2UV:  out Vec from gp);
+		     
+    	---Purpose: Returns all first and second derivatives.
+    
     
 
     Tolerance(myclass; Quad: Quadric from IntSurf )

src/IntSurf/IntSurf_QuadricTool.cxx

@@ -17,12 +17,47 @@
 #include <gp_Sphere.hxx>
 #include <gp_Cylinder.hxx>
 
+
+void IntSurf_QuadricTool::ValueAndGradient (const IntSurf_Quadric& Quad,
+                                            const Standard_Real X,
+                                            const Standard_Real Y,
+                                            const Standard_Real Z,
+                                            Standard_Real& Val,
+                                            gp_Vec& Grad,
+                                            gp_Vec& D1U,
+                                            gp_Vec& D1V)
+{
+  gp_Pnt thePoint(X,Y,Z);
+  Quad.ValAndGrad(thePoint, Val, Grad);
+  Standard_Real U, V;
+  Quad.Parameters(thePoint, U, V);
+  Quad.D1(U, V, thePoint, D1U, D1V);
+}
+
+void IntSurf_QuadricTool::D2 (const IntSurf_Quadric& Quad,
+                              const gp_Pnt& aPoint,
+                              gp_Vec& D1U,
+                              gp_Vec& D1V,
+                              gp_Vec& D2U,
+                              gp_Vec& D2V,
+                              gp_Vec& D2UV)
+{
+  Standard_Real U, V;
+  Quad.Parameters(aPoint, U, V);
+  gp_Pnt thePoint;
+  Quad.D1(U, V, thePoint, D1U, D1V);
+  D2U = Quad.DN(U, V, 2, 0);
+  D2V = Quad.DN(U, V, 0, 2);
+  D2UV = Quad.DN(U, V, 1, 1);
+}
+
 Standard_Real IntSurf_QuadricTool::Tolerance (const IntSurf_Quadric& Q) {
   switch (Q.TypeQuadric()) {
   case GeomAbs_Sphere:
     return 2.e-6*Q.Sphere().Radius();
   case GeomAbs_Cylinder:
-    return 2.e-6*Q.Cylinder().Radius();
+    //return 2.e-6*Q.Cylinder().Radius();
+    return Precision::Confusion();
   default: 
     break;
   }

src/IntSurf/IntSurf_QuadricTool.lxx

@@ -31,15 +31,3 @@ inline void IntSurf_QuadricTool::Gradient (const IntSurf_Quadric& Quad,
 
   V = Quad.Gradient(gp_Pnt(X,Y,Z));
 }
-
-inline void IntSurf_QuadricTool::ValueAndGradient (const IntSurf_Quadric& Quad,
-						   const Standard_Real X,
-						   const Standard_Real Y,
-						   const Standard_Real Z,
-						   Standard_Real& Val,
-						   gp_Vec& V) {
-
-  Quad.ValAndGrad(gp_Pnt(X,Y,Z),Val,V);
-}
-
-

src/IntWalk/IntWalk.cdl

@@ -35,7 +35,7 @@ is
 
     enumeration StatusDeflection is 
                 PasTropGrand, PointConfondu, ArretSurPointPrecedent,
-                ArretSurPoint, OK;
+                ArretSurPoint, OK, OKtangent;
 
 -- StepTooGreat, ConfusedPoint, StopOnPreviousPoint, StopOnPoint, OK
     

src/IntWalk/IntWalk_IWLine.cdl

@@ -124,6 +124,11 @@ is
     	---C++: inline
     	is static;
 
+    SetTangency(me: mutable; IsTangency: Boolean from Standard)
+    
+    	---C++: inline
+    	is static;
+
 
     NbPoints(me)
     
@@ -307,6 +312,13 @@ is
 	
 	is static;
 
+    IsTangency(me)
+    
+    	returns Boolean from Standard
+    	---C++: inline
+	
+	is static;
+
 
 fields
 
@@ -323,5 +335,6 @@ fields
     vcttg         : Vec              from gp;
     istgtbeg      : Boolean          from Standard;
     istgtend      : Boolean          from Standard;
+    istangency    : Boolean          from Standard;
   
 end IWLine;

src/IntWalk/IntWalk_IWLine.gxx

@@ -20,7 +20,8 @@ IntWalk_IWLine::IntWalk_IWLine (const IntSurf_Allocator& theAllocator) :
   hasFirst (Standard_False), hasLast (Standard_False),
   firstIndex (-1), lastIndex (-1),
   indextg (-1),
-  istgtbeg (Standard_False), istgtend (Standard_False)
+  istgtbeg (Standard_False), istgtend (Standard_False),
+  istangency (Standard_False)
 {
 }
 

src/IntWalk/IntWalk_IWLine.lxx

@@ -162,6 +162,12 @@ inline void IntWalk_IWLine::SetTangencyAtEnd
   istgtend = IsTangent;
 }
 
+inline void IntWalk_IWLine::SetTangency
+  (const Standard_Boolean IsTangency) {
+
+  istangency = IsTangency;
+}
+
 inline const gp_Vec& IntWalk_IWLine::TangentVector
   (Standard_Integer& Index) const {
 //--    if(istgtend == Standard_False && istgtbeg == Standard_False) { 
@@ -181,4 +187,8 @@ inline Standard_Boolean IntWalk_IWLine::IsTangentAtEnd () const {
   return istgtend;
 }
 
+inline Standard_Boolean IntWalk_IWLine::IsTangency () const {
+
+  return istangency;
+}
 

src/IntWalk/IntWalk_IWalking.cdl

@@ -88,12 +88,15 @@ is
     Perform(me: in out;
             Pnts1    : ThePOPIterator;
             Pnts2    : ThePOLIterator;
+	    Pnts3    : ThePOLIterator;
             Func     : in out TheIWFunction;
 	    S        : ThePSurface;
             Reversed : Boolean from Standard = Standard_False)
      
-    	---Purpose: Searches a set of polylines starting on a point of Pnts1
+    	---Purpose: Searches a set of polylines starting on a point of
-    	--          or Pnts2.
+    	--          Pnts1 (border start points),
+    	--          or Pnts2 (inner start points),
+    	--          or Pnts3 (inner tangent start points)
     	--          Each point on a resulting polyline verifies F(u,v)=0
 
         is static;
@@ -180,7 +183,6 @@ is
 
         is static;
     
-
 -- -- private
 
     Cadrage(me; BornInf, BornSup, UVap : in out Vector from math;
@@ -206,16 +208,25 @@ is
     returns Boolean from Standard
     is static protected;		      
 
+    TestArretPassageTang (me: in out; Umult  : SequenceOfReal from TColStd;
+                                  Vmult  : SequenceOfReal from TColStd;
+                                  UV     : Vector from math;
+                                  Index  : Integer from Standard;
+	          	          Irang  : out Integer from Standard) 
+    returns Boolean from Standard
+    is static protected;		      
+
     TestArretAjout(me: in out; Section: in out TheIWFunction;
 		               UV     : in out Vector from math;
     	                       Irang  : out Integer from Standard; 
-                               PSol   : out PntOn2S from IntSurf )
+                               PSol   : out PntOn2S from IntSurf;
+			       OnPrevTangency : out Boolean from Standard)
     returns Boolean from Standard
     is static protected;		      
 		    
     TestArretCadre(me : in out; Umult  : SequenceOfReal from TColStd;
                                 Vmult  : SequenceOfReal from TColStd;
-                                Line   : TheIWLine;
+                                Line   : mutable TheIWLine;
 				Section: in out TheIWFunction;
                                 UV     : in out Vector from math;
                                 Irang  : out Integer from Standard)
@@ -227,9 +238,12 @@ is
                                UV              : Vector from math;
 		               StatusPrecedent : StatusDeflection from IntWalk;
                                NbDivision      : in out Integer from Standard;
-                               Step           : in out Real from Standard;
+                               Step            : in out Real from Standard;
 --                               StepV           : in out Real from Standard;
-                               StepSign        : Integer from Standard)
+                               StepSign        : Integer from Standard;
+			       CurNbPoints     : Integer from Standard;
+			       ArretAjout      : Boolean from Standard;
+			       OnPrevTangency  : Boolean from Standard)
     returns StatusDeflection from IntWalk		    
     is static protected;		      
 
@@ -244,20 +258,21 @@ is
                           Psol   : PntOn2S from IntSurf;
 			  Pnts1  : ThePOPIterator;
  			  Section: in out TheIWFunction;
-                          Line   : TheIWLine)
+                          Line   : mutable TheIWLine)
     is static protected;	
 
     ComputeCloseLine(me : in out; Umult  : SequenceOfReal from TColStd;
                                   Vmult  : SequenceOfReal from TColStd;
                                   Pnts1  : ThePOPIterator;
 				  Pnts2  : ThePOLIterator;
+				  Pnts3  : ThePOLIterator;
                                   Section: in out TheIWFunction;
                                   Rajout : in out Boolean from Standard)
     is static protected;	
 
     AddPointInCurrentLine(me ; N           : Integer from Standard;
                                PathPnt     : ThePointOfPath;
-                               CurrentLine : TheIWLine)
+                               CurrentLine : mutable TheIWLine)
     is static protected;
 
     MakeWalkingPoint(me : in out ; Case    : Integer from Standard;
@@ -267,6 +282,26 @@ is
 
     is static protected;
     
+    FindExactCuspPoint(me : in out ; Section : in out TheIWFunction;
+    	    	    	    	     Psol    : in out PntOn2S from IntSurf;
+				     PrevSignOfBcoeff : Integer from Standard;
+				     SignOfBcoeff     : Integer from Standard)
+    is static protected;
+    
+    FindExactTangentPoint(me : in out ; TolTang : Real from Standard;
+    	    	    	    	    	Section : in out TheIWFunction;
+    	    	    	    	    	Psol    : in out PntOn2S from IntSurf)
+    is static protected;
+    
+    ComputeDirOfTangentialIntersection(me : in out; Section : in out TheIWFunction;
+    	    	    	    	     	    	    StepSign: in out Integer from Standard;
+						    IsDiscriminantNull: out Boolean from Standard;
+						    SignOfBcoeff: in out Integer from Standard;
+						    NewD3d: out Vec   from gp;
+						    NewD2d: out Dir2d from gp)
+    returns Boolean from Standard
+    is static protected;
+    
     Clear (me: in out) is static protected;
     ---Purpose: Clears up internal containers
 
@@ -283,6 +318,7 @@ fields
 
     wd1                 : VectorOfWalkingData from IntWalk;
     wd2                 : VectorOfWalkingData from IntWalk;
+    wd3                 : VectorOfWalkingData from IntWalk;
     nbMultiplicities    : VectorOfInteger   from IntWalk;
     Um                  : Real              from Standard; -- Min U de la surf
     UM                  : Real              from Standard; -- Max U de la surf
@@ -291,7 +327,9 @@ fields
     previousPoint       : PntOn2S           from IntSurf;
     previousd3d         : Vec               from gp;
     previousd2d         : Dir2d             from gp;
+    IsTangentialIntersection : Boolean      from Standard;
     seqAjout            : SequenceOfInteger from TColStd; 
     lines               : SequenceOfIWLine;
+    
 end IWalking;
 

src/IntWalk/IntWalk_IWalking_1.gxx

@@ -20,6 +20,7 @@ OSD_Chronometer Chronrsnld;
 
 #include <NCollection_IncAllocator.hxx>
 #include <Precision.hxx>
+#include <TColStd_Array1OfBoolean.hxx>
 
 IntWalk_IWalking::IntWalk_IWalking (const Standard_Real Epsilon,
                                     const Standard_Real Deflection,
@@ -31,7 +32,9 @@ IntWalk_IWalking::IntWalk_IWalking (const Standard_Real Epsilon,
       epsilon(Epsilon*Epsilon),
       wd1 (IntWalk_VectorOfWalkingData::allocator_type (new NCollection_IncAllocator)),
       wd2 (wd1.get_allocator()),
-      nbMultiplicities (wd1.get_allocator())
+      wd3 (wd2.get_allocator()),
+      nbMultiplicities (wd1.get_allocator()),
+      IsTangentialIntersection(Standard_False)
 {
 }
 
@@ -48,11 +51,13 @@ void IntWalk_IWalking::Clear()
 {
   wd1.clear();
   wd2.clear();
+  wd3.clear();
   IntWalk_WalkingData aDummy;
   aDummy.etat = -10;
   aDummy.ustart = aDummy.vstart = 0.;
   wd1.push_back (aDummy);
   wd2.push_back (aDummy);
+  wd3.push_back (aDummy);
   nbMultiplicities.clear();
   nbMultiplicities.push_back (-1);
   
@@ -80,6 +85,7 @@ void IntWalk_IWalking::Clear()
 //
 void IntWalk_IWalking::Perform(const ThePOPIterator& Pnts1,
 			       const ThePOLIterator& Pnts2,
+                               const ThePOLIterator& Pnts3,
 			       TheIWFunction& Func,
 			       const ThePSurface& Caro,
 			       const Standard_Boolean Reversed)
@@ -90,6 +96,7 @@ void IntWalk_IWalking::Perform(const ThePOPIterator& Pnts1,
   Standard_Boolean Rajout = Standard_False;
   Standard_Integer nbPnts1 = Pnts1.Length();
   Standard_Integer nbPnts2 = Pnts2.Length();
+  Standard_Integer nbPnts3 = Pnts3.Length();
   Standard_Real U,V;
 
   Clear();
@@ -137,6 +144,14 @@ void IntWalk_IWalking::Perform(const ThePOPIterator& Pnts1,
     wd2.push_back (aWD2);
   }
 
+  wd3.reserve (nbPnts3);
+  for (I = 1; I <= nbPnts3; I++) {
+    IntWalk_WalkingData aWD3;
+    aWD3.etat = 13;
+    ThePointOfLoopTool::Value2d(Pnts3.Value(I), aWD3.ustart, aWD3.vstart);
+    wd3.push_back (aWD3);
+  }
+
   tolerance(1) = ThePSurfaceTool::UResolution(Caro,Precision::Confusion());
   tolerance(2) = ThePSurfaceTool::VResolution(Caro,Precision::Confusion());
 
@@ -160,15 +175,137 @@ void IntWalk_IWalking::Perform(const ThePOPIterator& Pnts1,
 
   // calculation of all open lines   
   if (nbPnts1 != 0)
-    ComputeOpenLine(Umult,Vmult,Pnts1,Func,Rajout); 
+    ComputeOpenLine(Umult,Vmult,Pnts1,Func,Rajout);
-
   // calculation of all closed lines 
-  if (nbPnts2 != 0)
+  if (nbPnts2 != 0 || nbPnts3 != 0)
-    ComputeCloseLine(Umult,Vmult,Pnts1,Pnts2,Func,Rajout); 
+    ComputeCloseLine(Umult,Vmult,Pnts1,Pnts2,Pnts3,Func,Rajout); 
-  
   for (I = 1; I <= nbPnts1; I++) { 
     if (wd1[I].etat >0) seqSingle.Append(Pnts1(I));
   }
+
+  //Split non-tangency lines by extremities of tangency lines
+  IntSurf_SequenceOfPntOn2S EndsOfTangencies;
+  Standard_Integer i, j, k;
+  const Standard_Real TolSq = Precision::Confusion() * Precision::Confusion();
+  for (i = 1; i <= lines.Length(); i++)
+    if (lines(i)->IsTangency())
+    {
+      const IntSurf_PntOn2S& anEndPointOn2S = lines(i)->Value(1);
+      const gp_Pnt& anEndPnt = anEndPointOn2S.Value();
+      Standard_Boolean exists = Standard_False;
+      for (j = 1; j <= EndsOfTangencies.Length(); j++)
+      {
+        const gp_Pnt& anEndOfTang = EndsOfTangencies(j).Value();
+        if (anEndOfTang.SquareDistance(anEndPnt) <= TolSq)
+        {
+          exists = Standard_True;
+          break;
+        }
+      }
+      if (!exists)
+        EndsOfTangencies.Append(anEndPointOn2S);
+      if (!lines(i)->IsClosed())
+      {
+        const IntSurf_PntOn2S& anEndPointOn2S = lines(i)->Value(lines(i)->NbPoints());
+        const gp_Pnt& anEndPnt = anEndPointOn2S.Value();
+        exists = Standard_False;
+        for (j = 1; j <= EndsOfTangencies.Length(); j++)
+        {
+          const gp_Pnt& anEndOfTang = EndsOfTangencies(j).Value();
+          if (anEndOfTang.SquareDistance(anEndPnt) <= TolSq)
+          {
+            exists = Standard_True;
+            break;
+          }
+        }
+        if (!exists)
+          EndsOfTangencies.Append(anEndPointOn2S);
+      }
+    }
+
+  if (!EndsOfTangencies.IsEmpty())
+  {
+    TColStd_Array1OfBoolean SplitFlags(1, EndsOfTangencies.Length());
+    SplitFlags.Init(Standard_False);
+    
+    for (;;)
+    {
+      for (i = 1; i <= lines.Length(); i++)
+        if (!lines(i)->IsTangency())
+        {
+          Standard_Boolean IsSplit = Standard_False;
+          for (j = 1; j <= EndsOfTangencies.Length(); j++)
+          {
+            if (SplitFlags(j))
+              continue;
+            //gp_Pnt2d endoftang = EndsOfTangencies(j).ValueOnSurface(reversed);
+            gp_Pnt endoftang = EndsOfTangencies(j).Value();
+            for (k = 2; k < lines(i)->NbPoints()-1; k++)
+            {
+              const IntSurf_PntOn2S& P1 = lines(i)->Value(k);
+              const IntSurf_PntOn2S& P2 = lines(i)->Value(k+1);
+              //gp_Pnt2d pt1 = P1.ValueOnSurface(reversed);
+              //gp_Pnt2d pt2 = P2.ValueOnSurface(reversed);
+              const gp_Pnt& pt1 = P1.Value();
+              const gp_Pnt& pt2 = P2.Value();
+              //gp_Vec2d V1(endoftang, pt1);
+              //gp_Vec2d V2(endoftang, pt2);
+              gp_Vec V1(endoftang, pt1);
+              gp_Vec V2(endoftang, pt2);
+              Standard_Real ScalProd = V1 * V2;
+              if (ScalProd < 0.)
+              {
+                SplitFlags(j) = Standard_True;
+                Standard_Boolean ToRemoveFirstPnt = Standard_False, ToRemoveSecondPnt = Standard_False;
+                /*
+                Standard_Real Dist = pt1.Distance(pt2);
+                Standard_Real DistToFirst   = endoftang.Distance(pt1);
+                Standard_Real DistToSecond  = endoftang.Distance(pt2);
+                if (DistToFirst < Max(2.*Precision::Confusion(), 0.5*Dist))
+                  ToRemoveFirstPnt = Standard_True;
+                if (DistToSecond < Max(2.*Precision::Confusion(), 0.5*Dist))
+                  ToRemoveSecondPnt = Standard_True;
+                */
+                Standard_Real SqDist = pt1.SquareDistance(pt2);
+                Standard_Real SqDistToFirst   = endoftang.SquareDistance(pt1);
+                Standard_Real SqDistToSecond  = endoftang.SquareDistance(pt2);
+                if (SqDistToFirst < Max(4.*TolSq, 0.25*SqDist))
+                  ToRemoveFirstPnt = Standard_True;
+                if (SqDistToSecond < Max(4.*TolSq, 0.25*SqDist))
+                  ToRemoveSecondPnt = Standard_True;
+                Handle(IntWalk_TheIWLine) aNewLine = new IntWalk_TheIWLine (new NCollection_IncAllocator());
+                aNewLine->SetTangencyAtBegining(Standard_True);
+                aNewLine->AddStatusFirst(Standard_False, Standard_False);
+                aNewLine->AddPoint(EndsOfTangencies(j));
+                for (Standard_Integer kk = (ToRemoveSecondPnt)? k+2 : k+1;
+                     kk <= lines(i)->NbPoints(); kk++)
+                  aNewLine->AddPoint(lines(i)->Value(kk));
+                if (lines(i)->HasLastPoint())
+                  aNewLine->AddStatusLast(Standard_True, lines(i)->LastPointIndex(), lines(i)->LastPoint());
+                else
+                  aNewLine->AddStatusLast(Standard_False);
+                aNewLine->SetTangencyAtEnd(lines(i)->IsTangentAtEnd());
+                lines.Append(aNewLine);
+                lines(i)->Cut((ToRemoveFirstPnt)? k : k+1);
+                lines(i)->AddPoint(EndsOfTangencies(j));
+                lines(i)->AddStatusLast(Standard_False);
+                lines(i)->SetTangencyAtEnd(Standard_True);
+                IsSplit = Standard_True;
+                break;
+              }
+            }
+            if (IsSplit)
+              break;
+          }
+          if (IsSplit)
+            break;
+        }
+      if (i > lines.Length())
+        break;
+    }
+  }
+  ///////////////////////////////////////////////////////////
+  
   done = Standard_True;
 }
 

src/IntWalk/IntWalk_IWalking_2.gxx

@@ -214,43 +214,45 @@ Standard_Boolean IntWalk_IWalking::TestArretPassage
     previousPoint.ParametersOnS1(Up,Vp);
   }
 
-  for (size_t i = 1; i < wd2.size(); i++) { 
+  IntWalk_VectorOfWalkingData wdata [2] = {wd2, wd3};
-    if (wd2[i].etat > 0) { 
+  for (Standard_Integer ind = 0; ind < 2; ind++)
-      // debug jag 05.04.94
+    for (size_t i = 1; i < wdata[ind].size(); i++) { 
-
+      if (wdata[ind][i].etat > 0) { 
-//      if ((Up-wd2[i].ustart)*(UV(1)-wd2[i].ustart) +
+        // debug jag 05.04.94
-//	  (Vp-wd2[i].vstart)*(UV(2)-wd2[i].vstart) <= 0)
+        
-      Utest = wd2[i].ustart;
+        //      if ((Up-wd2[i].ustart)*(UV(1)-wd2[i].ustart) +
-      Vtest = wd2[i].vstart;
+        //	  (Vp-wd2[i].vstart)*(UV(2)-wd2[i].vstart) <= 0)
-
+        Utest = wdata[ind][i].ustart;
-      Du  = UV(1)-Utest;
+        Vtest = wdata[ind][i].vstart;
-      Dv  = UV(2)-Vtest;
+        
-      Dup = Up - Utest;
+        Du  = UV(1)-Utest;
-      Dvp = Vp - Vtest;
+        Dv  = UV(2)-Vtest;
-      
+        Dup = Up - Utest;
-//-- lbr le 30 oct 97 
+        Dvp = Vp - Vtest;
-
+        
-      //IFV for OCC20285
+        //-- lbr le 30 oct 97 
-
+        
-      if ((Abs(Du) < tolu2 && Abs(Dv) < tolv2) ||
+        //IFV for OCC20285
-	  (Abs(Dup) < tolu2 && Abs(Dvp) < tolv2)) { 
+        
-	    
+        if ((Abs(Du) < tolu2 && Abs(Dv) < tolv2) ||
-	wd2[i].etat = -wd2[i].etat;
+            (Abs(Dup) < tolu2 && Abs(Dvp) < tolv2)) { 
-      }
+          
-      else {
+          wdata[ind][i].etat = -wdata[ind][i].etat;
-	Standard_Real DDu = (UV(1)-Up);
+        }
-	Standard_Real DDv = (UV(2)-Vp);
+        else {
-	Standard_Real DDD = DDu*DDu+DDv*DDv;
+          Standard_Real DDu = (UV(1)-Up);
-	Standard_Real DD1 = Du*Du+Dv*Dv;
+          Standard_Real DDv = (UV(2)-Vp);
-	if(DD1<=DDD) { 
+          Standard_Real DDD = DDu*DDu+DDv*DDv;
-	  Standard_Real DD2 = Dup*Dup+Dvp*Dvp;
+          Standard_Real DD1 = Du*Du+Dv*Dv;
-	  if(DD2<=DDD && ((Du*Dup) + (Dv*Dvp*tolu/tolv) <= 0.)) { 	
+          if(DD1<=DDD) { 
-	    wd2[i].etat = -wd2[i].etat;
+            Standard_Real DD2 = Dup*Dup+Dvp*Dvp;
-	  }
+            if(DD2<=DDD && ((Du*Dup) + (Dv*Dvp*tolu/tolv) <= 0.)) { 	
-	}
+              wdata[ind][i].etat = -wdata[ind][i].etat;
+            }
+          }
+        }
       }
     }
-  }
 
   // stop test on point given at input and not yet processed
 
@@ -473,6 +475,62 @@ Standard_Boolean IntWalk_IWalking::TestArretPassage
       }
     }
   }
+  for (k = 1; k < (int)wd3.size(); k++) { 
+    if (wd3[k].etat > 0) {
+      Utest = wd3[k].ustart;
+      Vtest = wd3[k].vstart;
+      
+      Utest/=deltau;
+      Vtest/=deltav;
+      
+      Standard_Real UV1mUtest=UV1-Utest;
+      Standard_Real UV2mVtest=UV2-Vtest;
+      if(   (UV1mUtest<tolu2 && UV1mUtest>-tolu2)
+	 && (UV2mVtest<tolv2 && UV2mVtest>-tolv2)) {
+        
+	  wd3[k].etat = -wd3[k].etat; //-- mark the point as a crossing point 
+      }
+      else { 
+	Standard_Real UpmUtest = (Up-Utest);
+	Standard_Real VpmVtest = (Vp-Vtest);
+	Standard_Real dPreviousStart = (UpmUtest)*(UpmUtest)+(VpmVtest)*(VpmVtest);
+	Standard_Real dCurrentStart  = UV1mUtest * UV1mUtest + UV2mVtest * UV2mVtest;
+
+	Scal=(UpmUtest)*(UV1mUtest)+(VpmVtest)*(UV2mVtest);
+	if( (Abs(UpmUtest)<tolu && Abs(VpmVtest)<tolv)) {
+          
+	    wd3[k].etat = -wd3[k].etat;
+	}
+	else if(Scal<0 && (dPreviousStart+dCurrentStart < dPreviousCurrent)) {
+          
+	    wd3[k].etat = -wd3[k].etat; // mark the point as a crossing point 
+	}
+	else {
+	  if(dPreviousStart < dPreviousCurrent*0.25) {
+            
+	    wd3[k].etat = -wd3[k].etat; // mark the point as a crossing point 
+	    //-- cout<<"**** etat2 : ("<<k<<")"<<endl;
+	  }
+	  else { 
+	    if(dCurrentStart < dPreviousCurrent*0.25) {
+	      //-- cout<<"***** etat2 : ("<<k<<")"<<endl;
+	      wd3[k].etat = -wd3[k].etat; // mark the point as a crossing point 
+	    }
+	    else { 
+	      Standard_Real UMidUtest = 0.5*(UV1+Up)-Utest;
+	      Standard_Real VMidVtest = 0.5*(UV2+Vp)-Vtest;	    
+	      Standard_Real dMiddleStart =  UMidUtest*UMidUtest + VMidVtest*VMidVtest;
+
+	      if(dMiddleStart < dPreviousCurrent*0.5) { 
+		//-- cout<<"*********** etat2 : ("<<k<<")"<<endl;
+		wd3[k].etat = -wd3[k].etat; // mark the point as a crossing point 
+	      }
+	    }
+	  }
+	}
+      }
+    }
+  }
 
   // crossing test on crossing points.
   
@@ -507,12 +565,172 @@ Standard_Boolean IntWalk_IWalking::TestArretPassage
   return Arrive;
 }
 
+Standard_Boolean IntWalk_IWalking::TestArretPassageTang
+  (const TColStd_SequenceOfReal& Umult,
+   const TColStd_SequenceOfReal& Vmult,
+   const math_Vector& UV, 
+   const Standard_Integer Index, 
+   Standard_Integer& Irang)
+{
+// Umult, Vmult : table of stop (or crossing) points on border, here
+//          only crossing points are taken into account.
+// UV     : the current point.
+// Index  : index of the start point in uvstart2 of the current line
+//          (this is an interior point).
+// Irang  : at output : gives the index of the point passing in uvstart1 or 0.
+//          consider that there is risk to cross only one crossing point.
+
+// test of stop for a CLOSED intersection line.
+// 1) test of crossing on all interior points.
+// 2) test of crossing on all crossing points.
+
+  Standard_Real Up, Vp, Scal;
+  Standard_Boolean Arrive = Standard_False;
+  Standard_Integer N, k, i;
+  Standard_Real Utest,Vtest;
+  Standard_Real tolu = tolerance(1);
+  Standard_Real tolv = tolerance(2);
+
+  
+  // tests of stop and of crossing on all interior points.
+
+  if (!reversed) {
+    previousPoint.ParametersOnS2(Up,Vp);
+  }
+  else {
+    previousPoint.ParametersOnS1(Up,Vp);
+  }
+
+  Standard_Real UV1=UV(1);
+  Standard_Real UV2=UV(2);
+
+
+  //-- Put everything in one box 0 1  x 0 1 
+  //-- actually it is necessary to carry out tests in 3D
+
+  Standard_Real deltau=UM-Um;
+  Standard_Real deltav=VM-Vm;
+
+  Up/=deltau; UV1/=deltau; 
+  Vp/=deltav; UV2/=deltav;
+
+  tolu/=deltau;
+  tolv/=deltav;
+
+  Standard_Real tolu2=tolu+tolu;
+  Standard_Real tolv2=tolv+tolv;
+
+  
+  Standard_Real dPreviousCurrent = (Up-UV1)*(Up-UV1)+(Vp-UV2)*(Vp-UV2);
+  for (k = 1; k < (int)wd3.size(); k++) { 
+    if (wd3[k].etat > 0) {
+      Utest = wd3[k].ustart;
+      Vtest = wd3[k].vstart;
+      
+      Utest/=deltau;
+      Vtest/=deltav;
+      
+      Standard_Real UV1mUtest=UV1-Utest;
+      Standard_Real UV2mVtest=UV2-Vtest;
+      if(   (UV1mUtest<tolu2 && UV1mUtest>-tolu2)
+	 && (UV2mVtest<tolv2 && UV2mVtest>-tolv2)) { 
+	if(Index!=k) { 
+	  //-- cout<<"* etat2 : ("<<k<<")"<<endl;
+	  wd3[k].etat = -wd3[k].etat; //-- mark the point as a crossing point 
+	}
+	else {  //-- Index == k
+	  //-- cout<<"* Arrive"<<endl;
+	  Arrive=Standard_True;
+	}
+      }
+      else { 
+	Standard_Real UpmUtest = (Up-Utest);
+	Standard_Real VpmVtest = (Vp-Vtest);
+	Standard_Real dPreviousStart = (UpmUtest)*(UpmUtest)+(VpmVtest)*(VpmVtest);
+	Standard_Real dCurrentStart  = UV1mUtest * UV1mUtest + UV2mVtest * UV2mVtest;
+
+	Scal=(UpmUtest)*(UV1mUtest)+(VpmVtest)*(UV2mVtest);
+	if( (Abs(UpmUtest)<tolu && Abs(VpmVtest)<tolv)) { 
+	  if(Index != k ) { 
+	    //-- cout<<"** etat2 : ("<<k<<")"<<endl;
+	    wd3[k].etat = -wd3[k].etat;
+	  }
+	}
+	else if(Scal<0 && (dPreviousStart+dCurrentStart < dPreviousCurrent)) { 
+	  if (Index == k ) { // on a boucle.
+	    Arrive = Standard_True;
+	    //-- cout<<"** Arrive  : k="<<k<<endl;
+	  }
+	  else {
+	    //-- cout<<"*** etat2 : ("<<k<<")"<<endl;
+	    wd3[k].etat = -wd3[k].etat; // mark the point as a crossing point 
+	  }
+	}
+	else if(k!=Index) {
+	  if(dPreviousStart < dPreviousCurrent*0.25) { 
+	    wd3[k].etat = -wd3[k].etat; // mark the point as a crossing point 
+	    //-- cout<<"**** etat2 : ("<<k<<")"<<endl;
+	  }
+	  else { 
+	    if(dCurrentStart < dPreviousCurrent*0.25) {
+	      //-- cout<<"***** etat2 : ("<<k<<")"<<endl;
+	      wd3[k].etat = -wd3[k].etat; // mark the point as a crossing point 
+	    }
+	    else { 
+	      Standard_Real UMidUtest = 0.5*(UV1+Up)-Utest;
+	      Standard_Real VMidVtest = 0.5*(UV2+Vp)-Vtest;	    
+	      Standard_Real dMiddleStart =  UMidUtest*UMidUtest + VMidVtest*VMidVtest;
+
+	      if(dMiddleStart < dPreviousCurrent*0.5) { 
+		//-- cout<<"*********** etat2 : ("<<k<<")"<<endl;
+		wd3[k].etat = -wd3[k].etat; // mark the point as a crossing point 
+	      }
+	    }
+	  }
+	}
+      }
+    }
+  }
+
+  // crossing test on crossing points.
+  
+  Irang =0;
+  for (i = 1; i < (int)wd1.size(); i++) {
+    if (wd1[i].etat > 0 && wd1[i].etat < 11) { //test of crossing points 
+      Utest = wd1[i].ustart;
+      Vtest = wd1[i].vstart;
+      Utest/=deltau;
+      Vtest/=deltav;
+      
+      if (((Up-Utest) * (UV1-Utest) + (Vp-Vtest) * (UV2-Vtest) < 0) ||
+	  (Abs(UV1-Utest) < tolu &&  Abs(UV2-Vtest) < tolv)) 
+	Irang = i;
+      else if (nbMultiplicities[i] > 0) {
+	N=0;
+	for (k = 1; k < i; k++) N = N + nbMultiplicities[k];
+	for (Standard_Integer j = N + 1; j <= N + nbMultiplicities[i]; j++) {
+	  Standard_Real Umultj = Umult(j)/deltau;
+	  Standard_Real Vmultj = Vmult(j)/deltav;	  
+          if (((Up-Umultj)*(UV1-Umultj) +
+	       (Vp-Vmultj)*(UV2-Vmultj) < 0) ||
+	      (Abs(UV1-Umultj) < tolu &&
+	       Abs(UV2-Vmultj) < tolv)) {
+	    Irang=i;
+	    break;
+          }
+	}
+      }
+    }    
+  }
+  return Arrive;
+}
 
 Standard_Boolean IntWalk_IWalking::TestArretAjout
   (TheIWFunction& sp,
    math_Vector& UV, 
    Standard_Integer& Irang,
-   IntSurf_PntOn2S& Psol) 
+   IntSurf_PntOn2S& Psol,
+   Standard_Boolean& OnPrevTangency)
 
 // test of stop on added points 
 // these are the points on the natural biorder that were not given at input
@@ -559,7 +777,8 @@ Standard_Boolean IntWalk_IWalking::TestArretAjout
 	  (Abs(UV(1)-U1) < tolerance(1) &&  
 	   Abs(UV(2)-V1) < tolerance(2))) {
 //jag 940615	Irang= -Abs(Irang); 
-	Arrive = Standard_True; 
+	Arrive = Standard_True;
+        OnPrevTangency = Line->IsTangency();
 	UV(1) = U1;
 	UV(2) = V1;
 	Standard_Real abidF[1], abidD[1][2];

src/IntWalk/IntWalk_IWalking_3.gxx

@@ -66,10 +66,11 @@ void IntWalk_IWalking::ComputeOpenLine(const TColStd_SequenceOfReal& Umult,
   Standard_Integer I, N;
   Standard_Real aBornInf[2], aBornSup[2], aUVap[2];
   math_Vector BornInf(aBornInf,1,2), BornSup(aBornSup,1,2), UVap(aUVap,1,2);
-  Standard_Real PasC, PasCu, PasCv;
+  Standard_Real PasC, PasCu, PasCv, PasSav;
   Standard_Boolean Arrive; // shows if the line ends
   Standard_Boolean Cadre;  // shows if one is on border of the domain
   Standard_Boolean ArretAjout;  //shows if one is on added point
+  Standard_Boolean OnPrevTangency; //shows if one is on added point of tangent line
   IntSurf_PntOn2S Psol;
   Handle(IntWalk_TheIWLine)  CurrentLine;    // line under construction
   Standard_Boolean Tgtend;
@@ -80,6 +81,9 @@ void IntWalk_IWalking::ComputeOpenLine(const TColStd_SequenceOfReal& Umult,
   // number of divisions of step for each section
 
   Standard_Integer StepSign;
+  Standard_Integer StepSignTangent;
+  Standard_Integer SignOfBcoeff = 0;
+  Standard_Integer PrevSignOfBcoeff = SignOfBcoeff;
   
   ThePointOfPath PathPnt;
 
@@ -105,11 +109,16 @@ void IntWalk_IWalking::ComputeOpenLine(const TColStd_SequenceOfReal& Umult,
     // modified by NIZHNY-MKK  Fri Oct 27 12:33:37 2000.BEGIN
     if ((wd1[I].etat > 11) || ((wd1[I].etat < -11) && (movementdirectioninfo[I]!=0))) {
     // modified by NIZHNY-MKK  Fri Oct 27 12:33:43 2000.END
-      PathPnt = Pnts1.Value(I);     
+      PathPnt = Pnts1.Value(I);
+#ifdef DRAW
+      cout<<"PathPnt("<<I<<"):"<<endl;
+#endif
       CurrentLine = new IntWalk_TheIWLine (new NCollection_IncAllocator());
       CurrentLine->SetTangencyAtBegining(Standard_False);
       Tgtend = Standard_False;
       CurrentLine->AddStatusFirst(Standard_False, Standard_True, I, PathPnt);
+      IsTangentialIntersection = Standard_False;
+      SignOfBcoeff = PrevSignOfBcoeff = 0;
       UVap(1) = wd1[I].ustart;
       UVap(2) = wd1[I].vstart;
       MakeWalkingPoint(11, UVap(1), UVap(2), Func, previousPoint);
@@ -118,27 +127,27 @@ void IntWalk_IWalking::ComputeOpenLine(const TColStd_SequenceOfReal& Umult,
       CurrentLine->AddPoint(previousPoint);
       // modified by NIZHNY-MKK  Fri Oct 27 12:34:32 2000.BEGIN
       if(movementdirectioninfo[I] !=0) {
-        if(movementdirectioninfo[I] < 0) {
+	if(movementdirectioninfo[I] < 0) {
-          StepSign = -1;
+	  StepSign = -1;
-          CurrentLine->SetTangentVector(previousd3d.Reversed(),1);
+	  CurrentLine->SetTangentVector(previousd3d.Reversed(),1);
-        } else {
+	} else {
-          StepSign = 1; 
+	  StepSign = 1; 
-          CurrentLine->SetTangentVector(previousd3d,1);
+	  CurrentLine->SetTangentVector(previousd3d,1);
-        }
+	}
       } else {
-        Standard_Real tyutuyt=ThePointOfPathTool::Direction3d(PathPnt) * previousd3d;
+	Standard_Real tyutuyt=ThePointOfPathTool::Direction3d(PathPnt) * previousd3d;
-        if( tyutuyt < 0) {
+	if( tyutuyt < 0) {
-          StepSign = -1;
+	  StepSign = -1;
-          CurrentLine->SetTangentVector(previousd3d.Reversed(),1);
+	  CurrentLine->SetTangentVector(previousd3d.Reversed(),1);
-        }
+	}
-        else {
+	else {
-          StepSign = 1; 
+	  StepSign = 1; 
-          CurrentLine->SetTangentVector(previousd3d,1);
+	  CurrentLine->SetTangentVector(previousd3d,1);
-        }
+	}
       }
       // modified by NIZHNY-MKK  Fri Oct 27 12:34:37 2000.END
 
-      //  Modified by Sergey KHROMOV - Tue Nov 20 10:41:45 2001 Begin
+//  Modified by Sergey KHROMOV - Tue Nov 20 10:41:45 2001 Begin
       wd1[I].etat = - abs(wd1[I].etat);
       movementdirectioninfo[I] = (movementdirectioninfo[I]==0) ? StepSign : 0;
 //  Modified by Sergey KHROMOV - Tue Nov 20 10:41:56 2001 End
@@ -146,300 +155,357 @@ void IntWalk_IWalking::ComputeOpenLine(const TColStd_SequenceOfReal& Umult,
       Standard_Real d2dx = Abs(previousd2d.X()); 
       Standard_Real d2dy = Abs(previousd2d.Y()); 
       if (d2dx < tolerance(1)) {
-        PasC = pas * (VM-Vm)/d2dy;
+	PasC = pas * (VM-Vm)/d2dy;
       }
       else if (d2dy < tolerance(2)) {
-        PasC = pas * (UM-Um)/d2dx;
+	PasC = pas * (UM-Um)/d2dx;
       }
       else {
-        PasC = pas * Min((UM-Um)/d2dx,(VM-Vm)/d2dy);
+	PasC = pas * Min((UM-Um)/d2dx,(VM-Vm)/d2dy);
       }
 
+      PasSav = PasC;
+
       Arrive = Standard_False;
       ArretAjout = Standard_False;
+      OnPrevTangency = Standard_False;
       NbDivision = 0;
       StatusPrecedent = IntWalk_OK;
       // modified by NIZHNY-MKK  Fri Oct 27 12:39:37 2000
       Standard_Integer IndexOfPathPointDoNotCheck=0;
-      Standard_Integer aNbIter = 10;
+      //Standard_Integer aNbIter = 10;
+
       while (!Arrive) { //    as one of stop tests is not checked
-        Cadre = Cadrage(BornInf,BornSup,UVap,PasC,StepSign);
+
-        //  Border?
+	Cadre = Cadrage(BornInf,BornSup,UVap,PasC,StepSign);
+	//  Border?
 
 #ifdef CHRONO
-        Chronrsnld.Start();
+	Chronrsnld.Start();
 #endif
 
-        Rsnld.Perform(Func,UVap,BornInf,BornSup);
+	Rsnld.Perform(Func,UVap,BornInf,BornSup);
 
 #ifdef CHRONO
-        Chronrsnld.Stop();
+	Chronrsnld.Stop();
 #endif
 
-        if (Cadre) {
+	if (Cadre) {
-          BornInf(1) = Um; BornSup(1) = UM; BornInf(2) = Vm; BornSup(2) = VM;
+	  BornInf(1) = Um; BornSup(1) = UM; BornInf(2) = Vm; BornSup(2) = VM;
-        }
+	}
-        if (Rsnld.IsDone()) {
+	if (Rsnld.IsDone()) {
-          if (Abs(Func.Root()) > Func.Tolerance()) {
+	  if (Abs(Func.Root()) > Func.Tolerance()) {
-            PasC = PasC / 2.0;
+	    PasC = PasC / 2.0;
-            PasCu = Abs(PasC*previousd2d.X());
+	    PasCu = Abs(PasC*previousd2d.X());
-            PasCv = Abs(PasC*previousd2d.Y());
+	    PasCv = Abs(PasC*previousd2d.Y());
-            if (PasCu <= tolerance(1) && PasCv <= tolerance(2)) {
+	    if (PasCu <= tolerance(1) && PasCv <= tolerance(2)) {
-              if (CurrentLine->NbPoints() == 1) break;
+	      if (CurrentLine->NbPoints() == 1) break;
-              Arrive = Standard_True;
+	      Arrive = Standard_True;
-              CurrentLine->AddStatusLast(Standard_False);
+	      CurrentLine->AddStatusLast(Standard_False);
-              Tgtend = Standard_True; // check
+	      Tgtend = Standard_True; // check
               Rajout = Standard_True;
-              seqAjout.Append(lines.Length() + 1);
+	      seqAjout.Append(lines.Length() + 1);
-            }  
+	    }  
-          }
+	  }
-          else { // test stop
+	  else { // test stop
-            Rsnld.Root(UVap);
+	    Rsnld.Root(UVap);
-            Arrive = TestArretPassage(Umult, Vmult, Func, UVap, N);
+	    Arrive = TestArretPassage(Umult, Vmult, Func, UVap, N);
-            if (Arrive) {
+	    if (Arrive) {
-              Cadre = Standard_False;
+	      Cadre = Standard_False;
-              //in case if there is a frame and arrive at the same time
+	      //in case if there is a frame and arrive at the same time
-            }
+	    }
-            else {
+	    else {
-              if (Rajout) {
+	      if (Rajout) {
-                ArretAjout =TestArretAjout(Func, UVap, N, Psol);
+		ArretAjout =TestArretAjout(Func, UVap, N, Psol, OnPrevTangency);
-                if (ArretAjout) {
+		if (ArretAjout) {
-                  // jag 940615
+		  // jag 940615
-                  Tgtend = lines.Value(N)->IsTangentAtEnd();
+		  Tgtend = lines.Value(N)->IsTangentAtEnd();
-                  N = -N;
+		  N = -N;
-                }
+		}
-              }
+	      }
-              // modified by NIZHNY-MKK  Thu Nov  2 15:09:08 2000.BEGIN
+	      // modified by NIZHNY-MKK  Thu Nov  2 15:09:08 2000.BEGIN
-              if(!(Rajout && ArretAjout)) {
+	      if(!(Rajout && ArretAjout)) {
-                Standard_Real prevUp, prevVp;
+		Standard_Real prevUp, prevVp;
-                if (!reversed) {
+		if (!reversed) {
-                  previousPoint.ParametersOnS2(prevUp, prevVp);
+		  previousPoint.ParametersOnS2(prevUp, prevVp);
-                }
+		}
-                else {
+		else {
-                  previousPoint.ParametersOnS1(prevUp, prevVp);
+		  previousPoint.ParametersOnS1(prevUp, prevVp);
-                }
+		}
-                Arrive = TestPassedSolutionWithNegativeState(wd1, Umult, Vmult, prevUp, prevVp,
+		Arrive = TestPassedSolutionWithNegativeState(wd1, Umult, Vmult, prevUp, prevVp,
-                  nbMultiplicities, tolerance, Func, UVap, N);		
+							     nbMultiplicities, tolerance, Func, UVap, N);		
-                if(Arrive) {
+		if(Arrive) {
-                  Cadre = Standard_False;
+		  Cadre = Standard_False;
-                }
+		}
-              }
+	      }
-              // modified by NIZHNY-MKK  Thu Nov  2 15:09:13 2000.END
+	      // modified by NIZHNY-MKK  Thu Nov  2 15:09:13 2000.END
-              if (!ArretAjout && Cadre) {
+	      if (!ArretAjout && Cadre) {
-                if (CurrentLine->NbPoints() == 1) break; // cancel the line
+		if (CurrentLine->NbPoints() == 1) break; // cancel the line
-                TestArretCadre(Umult, Vmult, CurrentLine, Func, UVap, N);
+		TestArretCadre(Umult, Vmult, CurrentLine, Func, UVap, N);
-                //		if (N == 0) {
+//		if (N == 0) {
-                if (N <= 0) { // jag 941017
+		if (N <= 0) { // jag 941017
-                  MakeWalkingPoint(2, UVap(1), UVap(2), Func, Psol);
+		  MakeWalkingPoint(2, UVap(1), UVap(2), Func, Psol);
-                  Tgtend = Func.IsTangent();
+		  Tgtend = Func.IsTangent();
                   N = -N;
+		}
+	      }
+	    }
+	    Status = TestDeflection(Func, Arrive, UVap, StatusPrecedent,
+				    NbDivision,PasC,StepSign,
+                                    CurrentLine->NbPoints(),ArretAjout,OnPrevTangency);
+	    StatusPrecedent = Status;
+
+            if (Status == IntWalk_ArretSurPoint && Func.IsTangent())
+            {
+              PrevSignOfBcoeff = SignOfBcoeff;
+              gp_Vec   newd3d;
+              gp_Dir2d newd2d;
+              Standard_Boolean IsDiscriminantNull;
+              if (ComputeDirOfTangentialIntersection(Func, StepSignTangent,
+                                                     IsDiscriminantNull, SignOfBcoeff,
+                                                     newd3d, newd2d))
+              {
+                if (SignOfBcoeff == 0) //point of branching
+                  Status = IntWalk_ArretSurPoint;
+                else
+                {
+                  if (PrevSignOfBcoeff == 0 ||
+                      SignOfBcoeff == PrevSignOfBcoeff)
+                  {
+                    Status = IntWalk_OKtangent;
+                    MakeWalkingPoint(2, UVap(1), UVap(2), Func, previousPoint);
+                    previousd3d = newd3d;
+                    previousd2d = newd2d;
+                    CurrentLine->AddPoint(previousPoint);
+                  }
+                  else //change of sign of Bcoeff
+                  {
+                    MakeWalkingPoint(1, UVap(1), UVap(2), Func, Psol);
+                    FindExactCuspPoint(Func, Psol,
+                                       PrevSignOfBcoeff, SignOfBcoeff);
+                    Psol.ParametersOnSurface(reversed, UVap(1), UVap(2));
+                    Status = IntWalk_ArretSurPoint;
+                  }
                 }
               }
             }
-            Status = TestDeflection(Func, Arrive, UVap, StatusPrecedent,
+            
-              NbDivision,PasC,StepSign);
+	    if (Status == IntWalk_PasTropGrand) {
-            StatusPrecedent = Status;
+	      Arrive = Standard_False;
-            if (Status == IntWalk_PasTropGrand) {
+	      ArretAjout = Standard_False;
-              Arrive = Standard_False;
+	      Tgtend = Standard_False; // jag 940615
-              ArretAjout = Standard_False;
+	      if (!reversed) {
-              Tgtend = Standard_False; // jag 940615
+		previousPoint.ParametersOnS2(UVap(1), UVap(2));
-              if (!reversed) {
+	      }
-                previousPoint.ParametersOnS2(UVap(1), UVap(2));
+	      else {
-              }
+		previousPoint.ParametersOnS1(UVap(1), UVap(2));
-              else {
+	      }
-                previousPoint.ParametersOnS1(UVap(1), UVap(2));
+	    }
-              }
+	    else if (ArretAjout || Cadre) {
+	      Arrive = Standard_True;
+	      CurrentLine->AddStatusLast(Standard_False);
+	      if (Status != IntWalk_ArretSurPointPrecedent) {
+		CurrentLine->AddPoint(Psol);                      
+	      }
+	      if (Cadre && N==0) {
+		Rajout = Standard_True;
+		seqAjout.Append(lines.Length()+1);
+	      }
+	    }
+	    else if (Status == IntWalk_ArretSurPointPrecedent &&
+                     !Func.IsTangent()) //continue the line because no tangent point
+            {
+              Status = IntWalk_OK;
+              StatusPrecedent = IntWalk_OK;
+              PasC = PasSav;
+              MakeWalkingPoint(2, UVap(1), UVap(2), Func, previousPoint);
+	      previousd3d = Func.Direction3d();
+	      previousd2d = Func.Direction2d();
+	      CurrentLine->AddPoint(previousPoint);
             }
-            else if (ArretAjout || Cadre) {
+	    else if (Status == IntWalk_ArretSurPointPrecedent) {
-              Arrive = Standard_True;
+	      if (CurrentLine->NbPoints() == 1) { //cancel the line
-              CurrentLine->AddStatusLast(Standard_False);
+		Arrive = Standard_False;
-              if (Status != IntWalk_ArretSurPointPrecedent) {
+		break;
-                CurrentLine->AddPoint(Psol);                      
+	      }
-              }
+	      Arrive = Standard_True;
-              if (Cadre && N==0) {
+	      Rajout = Standard_True;
-                Rajout = Standard_True;
-                seqAjout.Append(lines.Length()+1);
-              }
-            }
-            else if (Status == IntWalk_ArretSurPointPrecedent) {
-              if (CurrentLine->NbPoints() == 1) { //cancel the line
-                Arrive = Standard_False;
-                break;
-              }
-              Arrive = Standard_True;
-              Rajout = Standard_True;
               seqAjout.Append(lines.Length() + 1);
-              CurrentLine->AddStatusLast(Standard_False);
+	      CurrentLine->AddStatusLast(Standard_False);
-              Tgtend = Standard_True; // check
+	      Tgtend = Standard_True; // check
-            }
+	    }
-            else if (Arrive)  {
+	    else if (Arrive)  {
-              if (CurrentLine->NbPoints() == 1 &&    // cancel the line
+	      if (CurrentLine->NbPoints() == 1 &&    // cancel the line
-                (N == I || Status == IntWalk_PointConfondu) ) {
+		  (N == I || Status == IntWalk_PointConfondu) ) {
-                  // if N == I the main uv is probably lost
+		// if N == I the main uv is probably lost
-                  // or the point is a point of accumulation
+		// or the point is a point of accumulation
-                  // if point is confused the start data is bad
+		// if point is confused the start data is bad
-                  Arrive =  Standard_False;
+		Arrive =  Standard_False;
-                  break;
+		break;
-              }
+	      }
-              // necessairily N > 0 jag 940617
+	      // necessairily N > 0 jag 940617
               // point of stop given at input 
-              PathPnt = Pnts1.Value(N);
+	      PathPnt = Pnts1.Value(N);
+	      
+	      Standard_Integer etat1N=wd1[N].etat;
+	      // modified by NIZHNY-MKK  Thu Nov  2 15:09:51 2000.BEGIN
+	      // 	      if (etat1N < 11) { // passing point that is a stop  
+	      if (Abs(etat1N) < 11) { // passing point that is a stop    
+		// modified by NIZHNY-MKK  Thu Nov  2 15:12:11 2000.END
+		if (Status == IntWalk_ArretSurPoint) { 
+		  CurrentLine->AddStatusLast(Standard_False);
+		  Tgtend = Standard_True; // need check
+		}
+		else { 
+		  Arrive = Standard_False;
+		}
+		CurrentLine->AddIndexPassing(N);
+	      }
+	      else { // point of stop given at input
+		if (etat1N == 11) {
+		  Tgtend = Standard_True;
+		}
+		CurrentLine->AddStatusLast(Standard_True, N, PathPnt);
+	      }
+	      AddPointInCurrentLine(N,PathPnt,CurrentLine);
+	      if ((etat1N != 1 && etat1N != 11)) {
+		// modified by NIZHNY-MKK  Fri Oct 27 12:43:05 2000.BEGIN
+		// 		wd1[N].etat= - wd1[N].etat;
+		wd1[N].etat = - Abs(etat1N);		
+		movementdirectioninfo[N] = (movementdirectioninfo[N]==0) ? StepSign : 0;
+		if(Arrive && movementdirectioninfo[N]!=0) {
+		  IndexOfPathPointDoNotCheck = N;
+		}
 
-              Standard_Integer etat1N=wd1[N].etat;
+		if(Arrive) {
-              // modified by NIZHNY-MKK  Thu Nov  2 15:09:51 2000.BEGIN
+		  Rajout = Standard_True;
-              // 	      if (etat1N < 11) { // passing point that is a stop  
+		  seqAjout.Append(lines.Length() + 1);
-              if (Abs(etat1N) < 11) { // passing point that is a stop    
+		}
-                // modified by NIZHNY-MKK  Thu Nov  2 15:12:11 2000.END
+		// modified by NIZHNY-MKK  Fri Oct 27 12:45:33 2000.END
-                if (Status == IntWalk_ArretSurPoint) { 
+	      }
-                  CurrentLine->AddStatusLast(Standard_False);
+	    }
-                  Tgtend = Standard_True; // need check
+	    else if (Status == IntWalk_ArretSurPoint) {
-                }
+              Arrive = Standard_True;
-                else { 
-                  Arrive = Standard_False;
-                }
-                CurrentLine->AddIndexPassing(N);
-              }
-              else { // point of stop given at input
-                if (etat1N == 11) {
-                  Tgtend = Standard_True;
-                }
-                CurrentLine->AddStatusLast(Standard_True, N, PathPnt);
-              }
-              AddPointInCurrentLine(N,PathPnt,CurrentLine);
-              if ((etat1N != 1 && etat1N != 11)) {
-                // modified by NIZHNY-MKK  Fri Oct 27 12:43:05 2000.BEGIN
-                // 		wd1[N].etat= - wd1[N].etat;
-                wd1[N].etat = - Abs(etat1N);		
-                movementdirectioninfo[N] = (movementdirectioninfo[N]==0) ? StepSign : 0;
-                if(Arrive && movementdirectioninfo[N]!=0) {
-                  IndexOfPathPointDoNotCheck = N;
-                }
-
-                if(Arrive) {
-                  Rajout = Standard_True;
-                  seqAjout.Append(lines.Length() + 1);
-                }
-                // modified by NIZHNY-MKK  Fri Oct 27 12:45:33 2000.END
-              }
-            }
-            else if (Status == IntWalk_ArretSurPoint) {
-              Arrive = Standard_True;                   
               CurrentLine->AddStatusLast(Standard_False);
               Tgtend = Standard_True;
               MakeWalkingPoint(1, UVap(1), UVap(2), Func, Psol);
               CurrentLine->AddPoint(Psol);
               Rajout = Standard_True;
               seqAjout.Append(lines.Length() + 1);
-            }
+	    }
-            else if (Status == IntWalk_OK) { 
+	    else if (Status == IntWalk_OK) { 
               MakeWalkingPoint(2, UVap(1), UVap(2), Func, previousPoint);
-              previousd3d = Func.Direction3d();
+	      previousd3d = Func.Direction3d();
-              previousd2d = Func.Direction2d();
+	      previousd2d = Func.Direction2d();
-              CurrentLine->AddPoint(previousPoint);
+	      CurrentLine->AddPoint(previousPoint);
-            }     
+	    }                           
             else if (Status == IntWalk_PointConfondu)
             {
-              aNbIter --;
+              //aNbIter --;
             }
-          }
+	  }
-        }
+	}
-        else { // no numerical solution
+	else { // no numerical solution
-          PasC = PasC / 2.;
+	  PasC = PasC / 2.;
-          PasCu = Abs(PasC*previousd2d.X());
+	  PasCu = Abs(PasC*previousd2d.X());
-          PasCv = Abs(PasC*previousd2d.Y());
+	  PasCv = Abs(PasC*previousd2d.Y());
-          if (PasCu <= tolerance(1) && PasCv <= tolerance(2)) {
+	  if (PasCu <= tolerance(1) && PasCv <= tolerance(2)) {
-            if (CurrentLine->NbPoints()==1) break;
+	    if (CurrentLine->NbPoints()==1) break;
-            Arrive = Standard_True;
+	    Arrive = Standard_True;
-            CurrentLine->AddStatusLast(Standard_False);
+	    CurrentLine->AddStatusLast(Standard_False);
-            Tgtend = Standard_True; // need check
+	    Tgtend = Standard_True; // need check
-            Rajout = Standard_True;
+	    Rajout = Standard_True;
             seqAjout.Append(lines.Length() + 1);
-          }  
+	  }  
-        }
+	}
 
+        /*
         if(aNbIter < 0)
           break;
+        */
       } // end of started line
       
       if (Arrive) {
-        CurrentLine->SetTangencyAtEnd(Tgtend);
+	CurrentLine->SetTangencyAtEnd(Tgtend);
-        lines.Append(CurrentLine);
+        CurrentLine->SetTangency(IsTangentialIntersection);
-        // modified by NIZHNY-MKK  Fri Oct 27 12:59:29 2000.BEGIN
+	lines.Append(CurrentLine);
-        movementdirectioninfo[I]=0;
+	// modified by NIZHNY-MKK  Fri Oct 27 12:59:29 2000.BEGIN
-        if(wd1[I].etat > 0)
+	movementdirectioninfo[I]=0;
-          // modified by NIZHNY-MKK  Fri Oct 27 12:59:42 2000.END
+	if(wd1[I].etat > 0)
-          wd1[I].etat=-wd1[I].etat;
+	// modified by NIZHNY-MKK  Fri Oct 27 12:59:42 2000.END
+	  wd1[I].etat=-wd1[I].etat;
 
-        //-- lbr le 5 juin 97 (Pb ds Contap)
+	//-- lbr le 5 juin 97 (Pb ds Contap)
-        for(Standard_Integer av=1; av<=nbPath; av++) { 
+	for(Standard_Integer av=1; av<=nbPath; av++) { 
-          // modified by NIZHNY-MKK  Fri Oct 27 13:00:22 2000.BEGIN
+	  // modified by NIZHNY-MKK  Fri Oct 27 13:00:22 2000.BEGIN
-          // 	  if (wd1[av].etat > 11) {
+	  // 	  if (wd1[av].etat > 11) {
-          if ((wd1[av].etat > 11) || 
+	  if ((wd1[av].etat > 11) || 
-            ((av!=I) && 
+	      ((av!=I) && 
-            (av!=IndexOfPathPointDoNotCheck) && 
+	       (av!=IndexOfPathPointDoNotCheck) && 
-            (wd1[av].etat < -11)  && 
+	       (wd1[av].etat < -11)  && 
-            (movementdirectioninfo[av]!=0)))
+	       (movementdirectioninfo[av]!=0))) {
-          {
+	  // modified by NIZHNY-MKK  Fri Oct 27 13:00:26 2000.END
-            // modified by NIZHNY-MKK  Fri Oct 27 13:00:26 2000.END
+	    Standard_Real Uav=wd1[av].ustart;
-            Standard_Real Uav=wd1[av].ustart;
+	    Standard_Real Vav=wd1[av].vstart;
-            Standard_Real Vav=wd1[av].vstart;
+	    Standard_Real Uavp,Vavp;
-            Standard_Real Uavp,Vavp;
+	    const IntSurf_PntOn2S &avP=CurrentLine->Value(CurrentLine->NbPoints());
-            const IntSurf_PntOn2S &avP=CurrentLine->Value(CurrentLine->NbPoints());
+	    if (!reversed) {
-            if (!reversed) {
+	      avP.ParametersOnS2(Uavp,Vavp);
-              avP.ParametersOnS2(Uavp,Vavp);
+	    }
-            }
+	    else {
-            else {
+	      avP.ParametersOnS1(Uavp,Vavp);
-              avP.ParametersOnS1(Uavp,Vavp);
+	    }
-            }
+	    Uav-=Uavp;
-            Uav-=Uavp;
+	    Vav-=Vavp;
-            Vav-=Vavp;
+	    Uav*=0.001; Vav*=0.001;
-            Uav*=0.001; Vav*=0.001;
+	    if(Abs(Uav)<tolerance(1) && Abs(Vav)<tolerance(2)) { 
-            if(Abs(Uav)<tolerance(1) && Abs(Vav)<tolerance(2)) { 
+	      // modified by NIZHNY-MKK  Fri Oct 27 13:01:38 2000.BEGIN
-              // modified by NIZHNY-MKK  Fri Oct 27 13:01:38 2000.BEGIN
+	      // 	      wd1[av].etat=-wd1[av].etat;
-              // 	      wd1[av].etat=-wd1[av].etat;
+	      if(wd1[av].etat < 0) {
-              if(wd1[av].etat < 0) {
+		movementdirectioninfo[av] = 0;
-                movementdirectioninfo[av] = 0;
+	      } else {
-              } else {
+		wd1[av].etat=-wd1[av].etat;
-                wd1[av].etat=-wd1[av].etat;
+		movementdirectioninfo[av] = StepSign;
-                movementdirectioninfo[av] = StepSign;
+	      }
-              }
+	      // modified by NIZHNY-MKK  Fri Oct 27 13:01:42 2000.END
-              // modified by NIZHNY-MKK  Fri Oct 27 13:01:42 2000.END
+	      CurrentLine->AddStatusLast(Standard_True, av, Pnts1.Value(av));
-              CurrentLine->AddStatusLast(Standard_True, av, Pnts1.Value(av));
+	      //-- cout<<"\n Debug ? lbr ds IntWalk_IWalking_3.gxx"<<endl;
-              //-- cout<<"\n Debug ? lbr ds IntWalk_IWalking_3.gxx"<<endl;
+	    }
-            }
+	    
-
+	    const IntSurf_PntOn2S &avPP=CurrentLine->Value(1);
-            const IntSurf_PntOn2S &avPP=CurrentLine->Value(1);
+	    if (!reversed) {
-            if (!reversed) {
+	      avPP.ParametersOnS2(Uavp,Vavp);
-              avPP.ParametersOnS2(Uavp,Vavp);
+	    }
-            }
+	    else {
-            else {
+	      avPP.ParametersOnS1(Uavp,Vavp);
-              avPP.ParametersOnS1(Uavp,Vavp);
+	    }
-            }
+	    Uav=wd1[av].ustart;
-            Uav=wd1[av].ustart;
+	    Vav=wd1[av].vstart;
-            Vav=wd1[av].vstart;
+	    Uav-=Uavp;
-            Uav-=Uavp;
+	    Vav-=Vavp;
-            Vav-=Vavp;
+	    Uav*=0.001; Vav*=0.001;
-            Uav*=0.001; Vav*=0.001;
+	    if(Abs(Uav)<tolerance(1) && Abs(Vav)<tolerance(2)) { 
-            if(Abs(Uav)<tolerance(1) && Abs(Vav)<tolerance(2)) { 
+	      // modified by NIZHNY-MKK  Fri Oct 27 13:02:49 2000.BEGIN
-              // modified by NIZHNY-MKK  Fri Oct 27 13:02:49 2000.BEGIN
+	      // 	      wd1[av].etat=-wd1[av].etat;
-              // 	      wd1[av].etat=-wd1[av].etat;
+	      if(wd1[av].etat < 0) {
-              if(wd1[av].etat < 0) {
+		movementdirectioninfo[av] = 0;
-                movementdirectioninfo[av] = 0;
+	      } else {
-              } else {
+		wd1[av].etat=-wd1[av].etat;
-                wd1[av].etat=-wd1[av].etat;
+		movementdirectioninfo[av] = -StepSign;
-                movementdirectioninfo[av] = -StepSign;
+	      }
-              }
+	      // modified by NIZHNY-MKK  Fri Oct 27 13:02:52 2000.END
-              // modified by NIZHNY-MKK  Fri Oct 27 13:02:52 2000.END
+	      //-- cout<<"\n Debug ? lbr ds IntWalk_IWalking_3.gxx"<<endl;
-              //-- cout<<"\n Debug ? lbr ds IntWalk_IWalking_3.gxx"<<endl;
+	      CurrentLine->AddStatusFirst(Standard_False, Standard_True, av, Pnts1.Value(av));
-              CurrentLine->AddStatusFirst(Standard_False, Standard_True, av, Pnts1.Value(av));
+	    }
-            }
+	  }
-          }
+	}
-        }
+	  
       }
     } //end of point processing
   } //end of all points
 }
 
+
 // modified by NIZHNY-MKK  Thu Nov  2 15:07:53 2000.BEGIN
 static Standard_Boolean TestPassedSolutionWithNegativeState(const IntWalk_VectorOfWalkingData& wd,
 							    const TColStd_SequenceOfReal& Umult,

src/IntWalk/IntWalk_IWalking_4.gxx

@@ -13,13 +13,17 @@
 // commercial license or contractual agreement.
 
 #include <NCollection_IncAllocator.hxx>
+#include <IntSurf_Quadric.hxx>
+#include <gp_Lin2d.hxx>
+#include <gce_MakeLin2d.hxx>
 
 void IntWalk_IWalking::ComputeCloseLine(const TColStd_SequenceOfReal& Umult,
 					const TColStd_SequenceOfReal& Vmult,
 					const ThePOPIterator& Pnts1,
 					const ThePOLIterator& Pnts2,
+					const ThePOLIterator& Pnts3,
 					TheIWFunction& Func,
-                                        Standard_Boolean& Rajout )
+                                        Standard_Boolean& Rajout ) 
 // *********** Processing of closed line **********************
 //
 // for any interior non-processed point 
@@ -41,7 +45,9 @@ void IntWalk_IWalking::ComputeCloseLine(const TColStd_SequenceOfReal& Umult,
 //
 // ******************************************************************** 
 {
-  Standard_Integer I,N = 0;
+  //const Standard_Real GlobalValTol = 1.e-12;
+
+  Standard_Integer I, N = 0;
   Standard_Real aBornInf[2], aBornSup[2], aUVap[2];
   math_Vector BornInf(aBornInf,1,2), BornSup(aBornSup,1,2);
   math_Vector Uvap(aUVap,1,2);// parameters of current approach
@@ -52,6 +58,7 @@ void IntWalk_IWalking::ComputeCloseLine(const TColStd_SequenceOfReal& Umult,
   Standard_Boolean Arrive;// show if line ends
   Standard_Boolean Cadre; // show if on border of the  domains
   Standard_Boolean ArretAjout; // show if on the added point
+  Standard_Boolean OnPrevTangency; //shows if one is on added point of tangent line
   IntSurf_PntOn2S Psol;
   Handle(IntWalk_TheIWLine)  CurrentLine; //line under construction
   ThePointOfPath PathPnt;
@@ -60,6 +67,9 @@ void IntWalk_IWalking::ComputeCloseLine(const TColStd_SequenceOfReal& Umult,
   Standard_Boolean Tgtbeg,Tgtend;
 
   Standard_Integer StepSign;
+  Standard_Integer StepSignTangent;
+  Standard_Integer SignOfBcoeff = 0;
+  Standard_Integer PrevSignOfBcoeff = SignOfBcoeff;
   
   IntWalk_StatusDeflection Status = IntWalk_OK, StatusPrecedent;
   Standard_Integer NbDivision ;   // number of divisions of step 
@@ -79,6 +89,7 @@ void IntWalk_IWalking::ComputeCloseLine(const TColStd_SequenceOfReal& Umult,
 
   // Check borders for degeneracy:
   math_Matrix D(1,1,1,2);
+  math_Vector FuncVal(1,1);
   const Standard_Integer aNbSamplePnt = 10;
   Standard_Boolean isLeftDegeneratedBorder[2] = {Standard_True, Standard_True};
   Standard_Boolean isRightDegeneratedBorder[2] = {Standard_True, Standard_True};
@@ -119,19 +130,30 @@ void IntWalk_IWalking::ComputeCloseLine(const TColStd_SequenceOfReal& Umult,
   }
 
   for (I = 1;I<=nbLoop;I++) {
-    if (wd2[I].etat > 12)
+    if (wd2[I].etat > 12) { // start point of closed line
-    { // start point of closed line
+      
       LoopPnt = Pnts2.Value(I);
-      previousPoint.SetValue(ThePointOfLoopTool::Value3d(LoopPnt),reversed,
+      gp_Pnt thePoint = ThePointOfLoopTool::Value3d(LoopPnt);
-           wd2[I].ustart,wd2[I].vstart);
+      previousPoint.SetValue(thePoint,reversed,
+			     wd2[I].ustart,wd2[I].vstart);
       previousd3d = ThePointOfLoopTool::Direction3d(LoopPnt);
       previousd2d = ThePointOfLoopTool::Direction2d(LoopPnt);
-
+      
+#ifdef DRAW
+      Standard_Real theU, theV;
+      ThePointOfLoopTool::Value2d(LoopPnt, theU, theV);
+      cout<<"LoopPnt("<<I<<"):"<<endl;
+      cout<<"u = "<<theU<<", v = "<<theV<<endl;
+      cout<<"pnt = "<<thePoint.X()<<" "<<thePoint.Y()<<" "<<thePoint.Z()<<endl<<endl;
+#endif
+      
       CurrentLine = new IntWalk_TheIWLine (new NCollection_IncAllocator());
       CurrentLine->AddPoint(previousPoint);
       CurrentLine->SetTangentVector(previousd3d,1);
       Tgtbeg = Standard_False;
       Tgtend = Standard_False;
+      IsTangentialIntersection = Standard_False;
+      SignOfBcoeff = PrevSignOfBcoeff = 0;
       Uvap(1) = wd2[I].ustart;
       Uvap(2) = wd2[I].vstart;
 
@@ -142,54 +164,54 @@ void IntWalk_IWalking::ComputeCloseLine(const TColStd_SequenceOfReal& Umult,
       Standard_Real d2dx = Abs(previousd2d.X()); 
       Standard_Real d2dy = Abs(previousd2d.Y()); 
       if (d2dx < tolerance(1)) {
-        PasC = pas * (VM-Vm)/d2dy;
+	PasC = pas * (VM-Vm)/d2dy;
       }
       else if (d2dy < tolerance(2)) {
-        PasC = pas * (UM-Um)/d2dx;
+	PasC = pas * (UM-Um)/d2dx;
       }
       else {
-        PasC = pas * Min((UM-Um)/d2dx,(VM-Vm)/d2dy);
+	PasC = pas * Min((UM-Um)/d2dx,(VM-Vm)/d2dy);
       }
 
       PasSav = PasC;
 
       Arrive = Standard_False;
       ArretAjout = Standard_False;
+      OnPrevTangency = Standard_False;
       NbDivision = 0;
       StatusPrecedent = IntWalk_OK;
       Standard_Integer aNbIter = 10;
       while (!Arrive) {  // as no test of stop is passed
-        Cadre=Cadrage(BornInf,BornSup,Uvap,PasC, StepSign);  // border?
+	Cadre=Cadrage(BornInf,BornSup,Uvap,PasC, StepSign);  // border?
 #ifdef CHRONO
-        Chronrsnld.Start();
+	Chronrsnld.Start();
 #endif
-
         Rsnld.Perform(Func,Uvap,BornInf,BornSup);
 
 #ifdef CHRONO
-        Chronrsnld.Stop();
+	Chronrsnld.Stop();
 #endif
 
-        if (Cadre) { // update of limits.
+	if (Cadre) { // update of limits.
-          BornInf(1) = Um;BornSup(1) = UM;BornInf(2) = Vm;BornSup(2) = VM;
+	  BornInf(1) = Um;BornSup(1) = UM;BornInf(2) = Vm;BornSup(2) = VM;
-        }
+	}
-        if (Rsnld.IsDone()) {
+	if (Rsnld.IsDone()) {
-          if (Abs(Func.Root()) > Func.Tolerance()) { // no solution for the tolerance
+	  if (Abs(Func.Root()) > Func.Tolerance()) { // no solution for the tolerance
-            PasC = PasC/2.;
+	    PasC = PasC/2.;
-            PasCu = Abs(PasC*previousd2d.X());
+	    PasCu = Abs(PasC*previousd2d.X());
-            PasCv = Abs(PasC*previousd2d.Y());
+	    PasCv = Abs(PasC*previousd2d.Y());
 
-            if (PasCu <= tolerance(1) && PasCv <= tolerance(2)) {
+	    if (PasCu <= tolerance(1) && PasCv <= tolerance(2)) {
-              if (CurrentLine->NbPoints()==1) break;
+	      if (CurrentLine->NbPoints()==1) break;
-              Arrive = Standard_True;
+	      Arrive = Standard_True;
-              CurrentLine->AddStatusFirstLast(Standard_False,
+	      CurrentLine->AddStatusFirstLast(Standard_False,
-                Standard_False,Standard_False);
+					      Standard_False,Standard_False);
-              Rajout = Standard_True;
+	      Rajout = Standard_True;
               seqAjout.Append(lines.Length()+1);
-              Tgtend = Standard_True;
+	      Tgtend = Standard_True;
-            }
+	    }
-          }
+	  }
-          else { // there is a solution
+	  else { // there is a solution
             Rsnld.Root(Uvap);
 
             // Avoid unitialized memory access.
@@ -231,199 +253,629 @@ void IntWalk_IWalking::ComputeCloseLine(const TColStd_SequenceOfReal& Umult,
               }
             }
 
-            Arrive = TestArretPassage(Umult,Vmult,Uvap,I,Ipass);
+	    Arrive = TestArretPassage(Umult,Vmult,Uvap,I,Ipass);
-            if (Arrive) {//reset proper parameter to test the arrow.
+	    if (Arrive) {//reset proper parameter to test the arrow.
-              Psol = CurrentLine->Value(1);
+	      Psol = CurrentLine->Value(1);
-              if (!reversed) {
+	      if (!reversed) {
-                Psol.ParametersOnS2(Uvap(1),Uvap(2));
+		Psol.ParametersOnS2(Uvap(1),Uvap(2));
-              }
+	      }
-              else {
+	      else {
-                Psol.ParametersOnS1(Uvap(1),Uvap(2));
+		Psol.ParametersOnS1(Uvap(1),Uvap(2));
-              }
+	      }
               Cadre=Standard_False; 
-              //in case if there is a frame and arrival at the same time
+	      //in case if there is a frame and arrival at the same time
+	    }
+	    else { // modif jag 940615
+
+	      if (Rajout) {    // test on added points
+		ArretAjout =TestArretAjout(Func, Uvap, N, Psol, OnPrevTangency);
+		if (ArretAjout) {
+		  if (N >0) {
+		    Tgtend = lines.Value(N)->IsTangentAtEnd();
+		    N = -N;
+		  }
+		  else {
+		    Tgtend = lines.Value(-N)->IsTangentAtBegining();
+		  }
+		  Arrive = (wd2[I].etat == 12);
+		}
+	      }
+
+	      if (!ArretAjout&& Cadre) {  // test on already marked points
+		if (CurrentLine->NbPoints() == 1)  break; // cancel the line
+		TestArretCadre(Umult,Vmult,CurrentLine,Func,Uvap,N);
+//		if (N==0) {
+		if (N <= 0) { // jag 941017
+		  MakeWalkingPoint(2,Uvap(1),Uvap(2),Func,Psol);
+		  Tgtend = Func.IsTangent(); // jag 940616
+		  N = -N;
+		}
+		Arrive = (wd2[I].etat == 12); // the line is open
+	      }
+	    }
+	    Status = TestDeflection(Func, Arrive,Uvap,StatusPrecedent,
+				    NbDivision,PasC,StepSign,
+                                    CurrentLine->NbPoints(),ArretAjout,OnPrevTangency);
+	    StatusPrecedent = Status; 
+	    if (Status == IntWalk_PasTropGrand) {// division of the step
+	      Arrive = Standard_False;
+	      ArretAjout = Standard_False;
+	      Tgtend = Standard_False; // jag 940616
+	      if (!reversed) {
+		previousPoint.ParametersOnS2(Uvap(1),Uvap(2));
+	      }
+	      else {
+		previousPoint.ParametersOnS1(Uvap(1),Uvap(2));
+	      }
+	    }
+	    else if (ArretAjout || Cadre) {
+
+	      if (Arrive) { // line s is open
+		CurrentLine->AddStatusLast(Standard_False);
+		if (Status != IntWalk_ArretSurPointPrecedent) {
+		  CurrentLine->AddPoint(Psol);                      
+		}
+		if (Cadre && N==0) {
+		  Rajout = Standard_True;
+		  seqAjout.Append(lines.Length()+1);
+		}
+                
+	      }
+	      else { // open
+		wd2[I].etat = 12; // declare it open
+		Tgtbeg = Tgtend;
+		Tgtend = Standard_False;
+		ArretAjout = Standard_False;
+		StepSign = -1;
+                StatusPrecedent = IntWalk_OK;
+		PasC = PasSav;
+		if (Status == IntWalk_ArretSurPointPrecedent) {
+		  OpenLine(0,Psol,Pnts1,Func,CurrentLine);
+		}
+		else {
+		  OpenLine(-lines.Length()-1,Psol,Pnts1,Func,CurrentLine);
+		}
+		if (Cadre && N==0) {
+		  Rajout = Standard_True;
+		  seqAjout.Append(-lines.Length()-1);
+		}
+	      }
+	    }
+	    else if (Status == IntWalk_ArretSurPointPrecedent &&
+                     !Func.IsTangent()) //continue the line because no tangent point
+            {
+              Status = IntWalk_OK;
+              StatusPrecedent = IntWalk_OK;
+              PasC = PasSav;
+
+	      if (Ipass!=0) CurrentLine->AddIndexPassing(Ipass);
+	      previousPoint.SetValue(Func.Point(),reversed,Uvap(1),Uvap(2));
+	      previousd3d = Func.Direction3d();
+	      previousd2d = Func.Direction2d();
+	      CurrentLine->AddPoint(previousPoint);
             }
-            else { // modif jag 940615
+	    else if ( Status == IntWalk_ArretSurPointPrecedent) {
-              if (Rajout) {    // test on added points
+	      if (CurrentLine->NbPoints() == 1) { //cancel the line
-                ArretAjout =TestArretAjout(Func,Uvap,N,Psol);
+		Arrive = Standard_False;
-                if (ArretAjout) {
+		break;
-                  if (N >0) {
+	      }
-                    Tgtend = lines.Value(N)->IsTangentAtEnd();
+	      if (wd2[I].etat >12) { //the line should become open
-                    N = -N;
+		wd2[I].etat = 12; //declare it open
+		ArretAjout = Standard_False;
+		OpenLine(0,Psol,Pnts1,Func,CurrentLine);
+		StepSign = -1;
+                StatusPrecedent = IntWalk_OK;
+		Arrive = Standard_False;
+		PasC = PasSav;
+		Rajout = Standard_True;
+                seqAjout.Append(-lines.Length()-1);
+	      }
+	      else { // line s is open                 
+		Arrive =Standard_True;
+		CurrentLine->AddStatusLast(Standard_False);
+		Rajout = Standard_True;
+                seqAjout.Append(lines.Length()+1);
+	      } 
+	    }
+	    else if (Arrive)  {
+	      if (wd2[I].etat > 12) {  //line closed good case
+		CurrentLine->AddStatusFirstLast(Standard_True,
+                                               Standard_False,Standard_False);
+		CurrentLine->AddPoint(CurrentLine->Value(1));              
+	      }
+	      else if (N >0) { //point of stop given at input 
+		PathPnt = Pnts1.Value(N);
+		CurrentLine->AddStatusLast(Standard_True,N,PathPnt);
+                AddPointInCurrentLine(N,PathPnt,CurrentLine);
+	      }
+	    }
+            else if (Status == IntWalk_ArretSurPoint)
+            {
+              if (Func.IsTangent())
+              {
+                PrevSignOfBcoeff = SignOfBcoeff;
+                gp_Vec   newd3d;
+                gp_Dir2d newd2d;
+                Standard_Boolean IsDiscriminantNull;
+                ComputeDirOfTangentialIntersection(Func, StepSignTangent,
+                                                   IsDiscriminantNull, SignOfBcoeff,
+                                                   newd3d, newd2d);
+                if (!IsDiscriminantNull) //we have gone from tangent line
+                {
+                  //take line from prevprevP to prevP
+                  //project current solution onto this line
+                  //direction of projection is direction to find null Discriminant again
+                  //and then repeat
+                  gp_Pnt2d prevprevP = CurrentLine->Value(CurrentLine->NbPoints()-1).ValueOnSurface(reversed);
+                  gp_Pnt2d prevP = previousPoint.ValueOnSurface(reversed);
+                  gp_Lin2d prevLine = gce_MakeLin2d(prevprevP, prevP);
+                  //project
+                  gp_Pnt2d CurSol(Uvap(1),Uvap(2));
+                  gp_XY LinLoc = prevLine.Location().XY();
+                  gp_XY LinDir = prevLine.Direction().XY();
+                  Standard_Real Parameter = (CurSol.XY() - LinLoc) * LinDir;
+                  gp_Pnt2d ProjCurSol( LinLoc + Parameter * LinDir );
+                  Uvap(1) = ProjCurSol.X();
+                  Uvap(2) = ProjCurSol.Y();
+                  Func.Values(Uvap, FuncVal, D);
+                  /////////
+                  ComputeDirOfTangentialIntersection(Func, StepSignTangent,
+                                                     IsDiscriminantNull, SignOfBcoeff,
+                                                     newd3d, newd2d);
+                }
+                if (SignOfBcoeff == 0) //point of branching
+                  Status = IntWalk_ArretSurPoint;
+                else
+                {
+                  if (PrevSignOfBcoeff == 0 ||
+                      SignOfBcoeff == PrevSignOfBcoeff)
+                  {                  
+                    Status = IntWalk_OKtangent;
+                    MakeWalkingPoint(2, Uvap(1), Uvap(2), Func, previousPoint);
+                    previousd3d = newd3d;
+                    previousd2d = newd2d;
+                    CurrentLine->AddPoint(previousPoint);
+                  }
+                  else //change of sign of Bcoeff
+                  {
+                    MakeWalkingPoint(1, Uvap(1), Uvap(2), Func, Psol);
+                    FindExactCuspPoint(Func, Psol,
+                                       PrevSignOfBcoeff, SignOfBcoeff);
+                    Psol.ParametersOnSurface(reversed, Uvap(1), Uvap(2));
+                    //Status = IntWalk_ArretSurPoint;
+                  }
+                }
+              }
+              if (Status == IntWalk_ArretSurPoint) {
+                if (wd2[I].etat >12) { //line should become open
+                  wd2[I].etat = 12; //declare it open
+                  Tgtbeg = Standard_True;
+                  Tgtend = Standard_False;
+                  N= -lines.Length()-1;
+                  Psol.SetValue(Func.Point(),reversed,Uvap(1),Uvap(2));
+                  OpenLine(N,Psol,Pnts1,Func,CurrentLine);
+                  StepSign = -1;
+                  SignOfBcoeff = PrevSignOfBcoeff = 0;
+                  Rajout = Standard_True;
+                  seqAjout.Append(N);
+                  StatusPrecedent = IntWalk_OK;
+                  Arrive = Standard_False;
+                  PasC = PasSav;	
+                }
+                else { 
+                  Arrive = Standard_True;                   
+                  if (Ipass!=0) { //point of passage, point of stop
+                    PathPnt = Pnts1.Value(Ipass);
+                    CurrentLine->AddStatusLast(Standard_True,Ipass,PathPnt);
+                    AddPointInCurrentLine(Ipass,PathPnt,CurrentLine);
                   }
                   else {
-                    Tgtend = lines.Value(-N)->IsTangentAtBegining();
+                    CurrentLine->AddStatusLast(Standard_False);
+                    IntSurf_PntOn2S newP;
+                    newP.SetValue(Func.Point(),reversed,Uvap(1),Uvap(2));
+                    CurrentLine->AddPoint(newP);
+                    Rajout = Standard_True;
+                    seqAjout.Append(lines.Length()+1);
                   }
-                  Arrive = (wd2[I].etat == 12);
-                }
-              }
-
-              if (!ArretAjout&& Cadre) {  // test on already marked points
-                if (CurrentLine->NbPoints() == 1)  break; // cancel the line
-                TestArretCadre(Umult,Vmult,CurrentLine,Func,Uvap,N);
-                //		if (N==0) {
-                if (N <= 0) { // jag 941017
-                  MakeWalkingPoint(2,Uvap(1),Uvap(2),Func,Psol);
-                  Tgtend = Func.IsTangent(); // jag 940616
-                  N = -N;
-                }
-                Arrive = (wd2[I].etat == 12); // the line is open
-              }
-            }
-            Status = TestDeflection(Func, Arrive,Uvap,StatusPrecedent,
-                                                NbDivision,PasC,StepSign);
-
-            StatusPrecedent = Status; 
-            if (Status == IntWalk_PasTropGrand) {// division of the step
-              Arrive = Standard_False;
-              ArretAjout = Standard_False;
-              Tgtend = Standard_False; // jag 940616
-              if (!reversed) {
-                previousPoint.ParametersOnS2(Uvap(1),Uvap(2));
-              }
-              else {
-                previousPoint.ParametersOnS1(Uvap(1),Uvap(2));
-              }
-            }
-            else if (ArretAjout || Cadre) {
-
-              if (Arrive) { // line s is open
-                CurrentLine->AddStatusLast(Standard_False);
-                if (Status != IntWalk_ArretSurPointPrecedent) {
-                  CurrentLine->AddPoint(Psol);                      
-                }
-                if (Cadre && N==0) {
-                  Rajout = Standard_True;
-                  seqAjout.Append(lines.Length()+1);
-                }
-
-              }
-              else { // open
-                wd2[I].etat = 12; // declare it open
-                Tgtbeg = Tgtend;
-                Tgtend = Standard_False;
-                ArretAjout = Standard_False;
-                StepSign = -1;
-                StatusPrecedent = IntWalk_OK;
-                PasC = PasSav;
-                if (Status == IntWalk_ArretSurPointPrecedent) {
-                  OpenLine(0,Psol,Pnts1,Func,CurrentLine);
-                }
-                else {
-                  OpenLine(-lines.Length()-1,Psol,Pnts1,Func,CurrentLine);
-                }
-                if (Cadre && N==0) {
-                  Rajout = Standard_True;
-                  seqAjout.Append(-lines.Length()-1);
                 }
               }
             }
-            else if ( Status == IntWalk_ArretSurPointPrecedent) {
+	    else if (Status == IntWalk_OK) { 
-              if (CurrentLine->NbPoints() == 1) { //cancel the line
+	      if (Ipass!=0) CurrentLine->AddIndexPassing(Ipass);
-                Arrive = Standard_False;
+	      previousPoint.SetValue(Func.Point(),reversed,Uvap(1),Uvap(2));
-                break;
+	      previousd3d = Func.Direction3d();
-              }
+	      previousd2d = Func.Direction2d();
-              if (wd2[I].etat >12) { //the line should become open
+	      CurrentLine->AddPoint(previousPoint);
-                wd2[I].etat = 12; //declare it open
+	    }
-                ArretAjout = Standard_False;
-                OpenLine(0,Psol,Pnts1,Func,CurrentLine);
-                StepSign = -1;
-                StatusPrecedent = IntWalk_OK;
-                Arrive = Standard_False;
-                PasC = PasSav;
-                Rajout = Standard_True;
-                seqAjout.Append(-lines.Length()-1);
-              }
-              else { // line s is open                 
-                Arrive =Standard_True;
-                CurrentLine->AddStatusLast(Standard_False);
-                Rajout = Standard_True;
-                seqAjout.Append(lines.Length()+1);
-              } 
-            }
-            else if (Arrive)  {
-              if (wd2[I].etat > 12) {  //line closed good case
-                CurrentLine->AddStatusFirstLast(Standard_True,
-                  Standard_False,Standard_False);
-                CurrentLine->AddPoint(CurrentLine->Value(1));              
-              }
-              else if (N >0) { //point of stop given at input 
-                PathPnt = Pnts1.Value(N);
-                CurrentLine->AddStatusLast(Standard_True,N,PathPnt);
-                AddPointInCurrentLine(N,PathPnt,CurrentLine);
-              }
-            }
-            else if (Status == IntWalk_ArretSurPoint) {
-              if (wd2[I].etat >12) { //line should become open
-                wd2[I].etat = 12; //declare it open
-                Tgtbeg = Standard_True;
-                Tgtend = Standard_False;
-                N= -lines.Length()-1;
-                Psol.SetValue(Func.Point(),reversed,Uvap(1),Uvap(2));
-                OpenLine(N,Psol,Pnts1,Func,CurrentLine);
-                StepSign = -1;
-                Rajout = Standard_True;
-                seqAjout.Append(N);
-                StatusPrecedent = IntWalk_OK;
-                Arrive = Standard_False;
-                PasC = PasSav;	
-              }
-              else { 
-                Arrive = Standard_True;                   
-                if (Ipass!=0) { //point of passage, point of stop
-                  PathPnt = Pnts1.Value(Ipass);
-                  CurrentLine->AddStatusLast(Standard_True,Ipass,PathPnt);
-                  AddPointInCurrentLine(Ipass,PathPnt,CurrentLine);
-                }
-                else {
-                  CurrentLine->AddStatusLast(Standard_False);
-                  IntSurf_PntOn2S newP;
-                  newP.SetValue(Func.Point(),reversed,Uvap(1),Uvap(2));
-                  CurrentLine->AddPoint(newP);
-                  Rajout = Standard_True;
-                  seqAjout.Append(lines.Length()+1);
-                }
-              }
-            }
-            else if (Status == IntWalk_OK) { 
-              if (Ipass!=0) CurrentLine->AddIndexPassing(Ipass);
-              previousPoint.SetValue(Func.Point(),reversed,Uvap(1),Uvap(2));
-              previousd3d = Func.Direction3d();
-              previousd2d = Func.Direction2d();
-              CurrentLine->AddPoint(previousPoint);
-            }
             else if (Status == IntWalk_PointConfondu)
             {
               aNbIter --;
             }
-          }
+	  }
-        }
+	}
-        else { //no numerical solution NotDone
+	else { //no numerical solution NotDone
-          PasC = PasC/2.;
+	  PasC = PasC/2.;
-          PasCu = Abs(PasC*previousd2d.X());
+	  PasCu = Abs(PasC*previousd2d.X());
-          PasCv = Abs(PasC*previousd2d.Y());
+	  PasCv = Abs(PasC*previousd2d.Y());
 
-          if (PasCu <= tolerance(1) && PasCv <= tolerance(2)) {
+	  if (PasCu <= tolerance(1) && PasCv <= tolerance(2)) {
             if (CurrentLine->NbPoints() == 1)  break; // cancel the line
-            Arrive = Standard_True;
+	    Arrive = Standard_True;
-            CurrentLine->AddStatusFirstLast(Standard_False,Standard_False,
+	    CurrentLine->AddStatusFirstLast(Standard_False,Standard_False,
-              Standard_False);
+                                           Standard_False);
-            Tgtend = Standard_True;
+	    Tgtend = Standard_True;
-            Rajout = Standard_True;
+	    Rajout = Standard_True;
-            seqAjout.Append(lines.Length()+1);
+	    seqAjout.Append(lines.Length()+1);
-          }  
+	  }  
-        }
+	}
-
+        
         if(aNbIter < 0)
           break;
       }// end of started line 
       if (Arrive) {
-        CurrentLine->SetTangencyAtBegining(Tgtbeg);
+	CurrentLine->SetTangencyAtBegining(Tgtbeg);
-        CurrentLine->SetTangencyAtEnd(Tgtend);
+	CurrentLine->SetTangencyAtEnd(Tgtend);
-
+        CurrentLine->SetTangency(IsTangentialIntersection);
-        lines.Append(CurrentLine);
+	
-        wd2[I].etat=-wd2[I].etat; //mark point as processed
+	lines.Append(CurrentLine);
+	wd2[I].etat=-wd2[I].etat; //mark point as processed
       }
     } //end of processing of start point
   } //end of all start points
+
+  //Process inner tangent points
+  IsTangentialIntersection = Standard_True;
+  StepSign = 1;
+  StepSignTangent = 1;
+  SignOfBcoeff = 0;
+  PrevSignOfBcoeff = SignOfBcoeff;
+  Standard_Integer nbTang = Pnts3.Length();
+  for (I = 1; I <= nbTang; I++) {
+    if (wd3[I].etat > 12) { // start point of closed line
+      
+      LoopPnt = Pnts3.Value(I);
+      gp_Pnt thePoint = ThePointOfLoopTool::Value3d(LoopPnt);
+      previousPoint.SetValue(thePoint,reversed,
+			     wd3[I].ustart, wd3[I].vstart);
+      Standard_Real theU, theV;
+      ThePointOfLoopTool::Value2d(LoopPnt, theU, theV);
+      IntSurf_PntOn2S PrevPointFromFunc;
+      MakeWalkingPoint(11, theU, theV, Func, PrevPointFromFunc);
+      /////
+      gp_Vec   newd3d;
+      gp_Dir2d newd2d;
+      Standard_Boolean IsDiscriminantNull;
+      if (!ComputeDirOfTangentialIntersection(Func, StepSignTangent,
+                                              IsDiscriminantNull, SignOfBcoeff,
+                                              newd3d, newd2d))
+        continue;
+      previousd3d = newd3d;
+      previousd2d = newd2d;
+      //here <previousd3d> and <previousd2d> are defined
+      
+      CurrentLine = new IntWalk_TheIWLine (new NCollection_IncAllocator());
+      CurrentLine->AddPoint(previousPoint);
+      CurrentLine->SetTangentVector(previousd3d,1);
+      Tgtbeg = Standard_False;
+      Tgtend = Standard_False;
+      SignOfBcoeff = PrevSignOfBcoeff = 0;
+      Uvap(1) = wd3[I].ustart;
+      Uvap(2) = wd3[I].vstart;
+
+      StepSign = 1;
+
+      // first step of advancement
+
+      Standard_Real d2dx = Abs(previousd2d.X()); 
+      Standard_Real d2dy = Abs(previousd2d.Y()); 
+      if (d2dx < tolerance(1)) {
+	PasC = pas * (VM-Vm)/d2dy;
+      }
+      else if (d2dy < tolerance(2)) {
+	PasC = pas * (UM-Um)/d2dx;
+      }
+      else {
+	PasC = pas * Min((UM-Um)/d2dx,(VM-Vm)/d2dy);
+      }
+
+      PasSav = PasC;
+
+      Arrive = Standard_False;
+      ArretAjout = Standard_False;
+      OnPrevTangency = Standard_False;
+      NbDivision = 0;
+      StatusPrecedent = IntWalk_OK;
+      while (!Arrive) {  // as no test of stop is passed
+	Cadre=Cadrage(BornInf,BornSup,Uvap,PasC, StepSign);  // border?
+#ifdef CHRONO
+	Chronrsnld.Start();
+#endif
+
+        //Standard_Boolean ExactStartSolution = Standard_False;
+        /*
+        for (Standard_Integer ind = 1; ind <= 2; ind++)
+        {
+          if (Uvap(ind) <= BornInf(ind)) Uvap(ind) = BornInf(ind);
+          else if (Uvap(ind) > BornSup(ind)) Uvap(ind) = BornSup(ind);
+        }
+        Func.Values(Uvap, FuncVal, D);
+        */
+        //if (FuncVal(1) >= GlobalValTol)
+        Rsnld.Perform(Func,Uvap,BornInf,BornSup);
+
+#ifdef CHRONO
+	Chronrsnld.Stop();
+#endif
+
+	if (Cadre) { // update of limits.
+	  BornInf(1) = Um;BornSup(1) = UM;BornInf(2) = Vm;BornSup(2) = VM;
+	}
+	if (Rsnld.IsDone()) {
+	  if (Abs(Func.Root()) > Func.Tolerance()) { // no solution for the tolerance
+	    PasC = PasC/2.;
+	    PasCu = Abs(PasC*previousd2d.X());
+	    PasCv = Abs(PasC*previousd2d.Y());
+
+	    if (PasCu <= tolerance(1) && PasCv <= tolerance(2)) {
+	      if (CurrentLine->NbPoints()==1) break;
+	      Arrive = Standard_True;
+	      CurrentLine->AddStatusFirstLast(Standard_False,
+					      Standard_False,Standard_False);
+	      Rajout = Standard_True;
+              seqAjout.Append(lines.Length()+1);
+	      Tgtend = Standard_True;
+	    }
+	  }
+	  else { // there is a solution
+            math_Vector NewUv(1,2), OldUv(1,2);
+            Rsnld.Root(NewUv);
+            gp_Pnt2d prevPnt2d, NewPnt2d(NewUv(1),NewUv(2)), Pnt2dap(Uvap(1),Uvap(2));
+            prevPnt2d = previousPoint.ValueOnSurface(reversed);
+            gp_Vec2d PrevToAp(prevPnt2d, Pnt2dap);
+            gp_Vec2d PrevToNew(prevPnt2d, NewPnt2d);
+            if (PrevToAp * PrevToNew < 0. ||
+                PrevToNew.Magnitude()/2. > PrevToAp.Magnitude()) //solution is bad
+            {
+              for (Standard_Integer ind = 1; ind <= 2; ind++)
+              {
+                if (Uvap(ind) <= BornInf(ind)) Uvap(ind) = BornInf(ind);
+                else if (Uvap(ind) > BornSup(ind)) Uvap(ind) = BornSup(ind);
+              }
+              Func.Values(Uvap, FuncVal, D);
+            }
+            else
+              Uvap = NewUv;
+              
+	    Arrive = TestArretPassageTang(Umult,Vmult,Uvap,I,Ipass);
+	    if (Arrive) {//reset proper parameter to test the arrow.
+	      Psol = CurrentLine->Value(1);
+	      if (!reversed) {
+		Psol.ParametersOnS2(Uvap(1),Uvap(2));
+	      }
+	      else {
+		Psol.ParametersOnS1(Uvap(1),Uvap(2));
+	      }
+              Cadre = Standard_False; 
+	      //in case if there is a frame and arrival at the same time
+	    }
+	    else { // modif jag 940615
+
+	      if (Rajout) {    // test on added points
+		ArretAjout = TestArretAjout(Func, Uvap, N, Psol, OnPrevTangency);
+		if (ArretAjout) {
+		  if (N >0) {
+		    Tgtend = lines.Value(N)->IsTangentAtEnd();
+		    N = -N;
+		  }
+		  else {
+		    Tgtend = lines.Value(-N)->IsTangentAtBegining();
+		  }
+		  Arrive = (wd3[I].etat == 12);
+		}
+	      }
+
+	      if (!ArretAjout&& Cadre) {  // test on already marked points
+		if (CurrentLine->NbPoints() == 1)  break; // cancel the line
+		TestArretCadre(Umult,Vmult,CurrentLine,Func,Uvap,N);
+//		if (N==0) {
+		if (N <= 0) { // jag 941017
+		  MakeWalkingPoint(2,Uvap(1),Uvap(2),Func,Psol);
+		  Tgtend = Func.IsTangent(); // jag 940616
+		  N = -N;
+		}
+		Arrive = (wd3[I].etat == 12); // the line is open
+	      }
+	    }
+	    Status = TestDeflection(Func, Arrive,Uvap,StatusPrecedent,
+				    NbDivision,PasC,StepSign,
+                                    CurrentLine->NbPoints(),ArretAjout,OnPrevTangency);
+	    StatusPrecedent = Status; 
+	    if (Status == IntWalk_PasTropGrand) {// division of the step
+	      Arrive = Standard_False;
+	      ArretAjout = Standard_False;
+	      Tgtend = Standard_False; // jag 940616
+	      if (!reversed) {
+		previousPoint.ParametersOnS2(Uvap(1),Uvap(2));
+	      }
+	      else {
+		previousPoint.ParametersOnS1(Uvap(1),Uvap(2));
+	      }
+	    }
+	    else if (ArretAjout || Cadre) {
+
+	      if (Arrive) { // line s is open
+		CurrentLine->AddStatusLast(Standard_False);
+		if (Status != IntWalk_ArretSurPointPrecedent) {
+		  CurrentLine->AddPoint(Psol);                      
+		}
+		if (Cadre && N==0) {
+		  Rajout = Standard_True;
+		  seqAjout.Append(lines.Length()+1);
+		}
+                
+	      }
+	      else { // open
+		wd3[I].etat = 12; // declare it open
+		Tgtbeg = Tgtend;
+		Tgtend = Standard_False;
+		ArretAjout = Standard_False;
+		StepSign = -1;
+                StatusPrecedent = IntWalk_OK;
+		PasC = PasSav;
+		if (Status == IntWalk_ArretSurPointPrecedent) {
+		  OpenLine(0,Psol,Pnts1,Func,CurrentLine);
+		}
+		else {
+		  OpenLine(-lines.Length()-1,Psol,Pnts1,Func,CurrentLine);
+		}
+		if (Cadre && N==0) {
+		  Rajout = Standard_True;
+		  seqAjout.Append(-lines.Length()-1);
+		}
+	      }
+	    }
+
+	    else if ( Status == IntWalk_ArretSurPointPrecedent) {
+	      if (CurrentLine->NbPoints() == 1) { //cancel the line
+		Arrive = Standard_False;
+		break;
+	      }
+	      if (wd3[I].etat >12) { //the line should become open
+		wd3[I].etat = 12; //declare it open
+		ArretAjout = Standard_False;
+		OpenLine(0,Psol,Pnts1,Func,CurrentLine);
+		StepSign = -1;
+                StatusPrecedent = IntWalk_OK;
+		Arrive = Standard_False;
+		PasC = PasSav;
+		Rajout = Standard_True;
+                seqAjout.Append(-lines.Length()-1);
+	      }
+	      else { // line s is open                 
+		Arrive = Standard_True;
+		CurrentLine->AddStatusLast(Standard_False);
+		Rajout = Standard_True;
+                seqAjout.Append(lines.Length()+1);
+	      } 
+	    }
+	    else if (Arrive)  {
+	      if (wd3[I].etat > 12) {  //line closed good case
+		CurrentLine->AddStatusFirstLast(Standard_True,
+                                               Standard_False,Standard_False);
+		CurrentLine->AddPoint(CurrentLine->Value(1));              
+	      }
+	      else if (N >0) { //point of stop given at input 
+		PathPnt = Pnts1.Value(N);
+		CurrentLine->AddStatusLast(Standard_True,N,PathPnt);
+                AddPointInCurrentLine(N,PathPnt,CurrentLine);
+	      }
+	    }
+            else if (Status == IntWalk_ArretSurPoint && Func.IsTangent())
+            {
+              PrevSignOfBcoeff = SignOfBcoeff;
+              gp_Vec   newd3d;
+              gp_Dir2d newd2d;
+              Standard_Boolean IsDiscriminantNull;
+              ComputeDirOfTangentialIntersection(Func, StepSignTangent,
+                                                 IsDiscriminantNull, SignOfBcoeff,
+                                                 newd3d, newd2d);
+              if (SignOfBcoeff == 0) //point of branching
+              {
+                //Status = IntWalk_ArretSurPoint;
+              }
+              else
+              {
+                if (PrevSignOfBcoeff == 0 ||
+                    SignOfBcoeff == PrevSignOfBcoeff)
+                {
+                  Status = IntWalk_OKtangent;
+                  MakeWalkingPoint(2, Uvap(1), Uvap(2), Func, previousPoint);
+                  previousd3d = newd3d;
+                  previousd2d = newd2d;
+                  CurrentLine->AddPoint(previousPoint);
+                }
+                else //change of sign of Bcoeff
+                {
+                  MakeWalkingPoint(1, Uvap(1), Uvap(2), Func, Psol);
+                  FindExactCuspPoint(Func, Psol,
+                                     PrevSignOfBcoeff, SignOfBcoeff);
+                  Psol.ParametersOnSurface(reversed, Uvap(1), Uvap(2));
+                  //Status = IntWalk_ArretSurPoint;
+                }
+              }
+              if (Status == IntWalk_ArretSurPoint) {
+                if (wd3[I].etat >12) { //line should become open
+                  wd3[I].etat = 12; //declare it open
+                  Tgtbeg = Standard_True;
+                  Tgtend = Standard_False;
+                  N= -lines.Length()-1;
+                  Psol.SetValue(Func.Point(),reversed,Uvap(1),Uvap(2));
+                  OpenLine(N,Psol,Pnts1,Func,CurrentLine);
+                  StepSign = -1;
+                  SignOfBcoeff = PrevSignOfBcoeff = 0;
+                  Rajout = Standard_True;
+                  seqAjout.Append(N);
+                  StatusPrecedent = IntWalk_OK;
+                  Arrive = Standard_False;
+                  PasC = PasSav;	
+                }
+                else { 
+                  Arrive = Standard_True;                   
+                  if (Ipass!=0) { //point of passage, point of stop
+                    PathPnt = Pnts1.Value(Ipass);
+                    CurrentLine->AddStatusLast(Standard_True,Ipass,PathPnt);
+                    AddPointInCurrentLine(Ipass,PathPnt,CurrentLine);
+                  }
+                  else {
+                    CurrentLine->AddStatusLast(Standard_False);
+                    IntSurf_PntOn2S newP;
+                    newP.SetValue(Func.Point(),reversed,Uvap(1),Uvap(2));
+                    CurrentLine->AddPoint(newP);
+                    Rajout = Standard_True;
+                    seqAjout.Append(lines.Length()+1);
+                  }
+                }
+              }
+            }
+	    else if (Status == IntWalk_OK) { 
+	      if (Ipass!=0) CurrentLine->AddIndexPassing(Ipass);
+	      previousPoint.SetValue(Func.Point(),reversed,Uvap(1),Uvap(2));
+	      //previousd3d = Func.Direction3d();
+	      //previousd2d = Func.Direction2d();
+	      CurrentLine->AddPoint(previousPoint);
+	    }
+	  }
+	}
+	else { //no numerical solution NotDone
+	  PasC = PasC/2.;
+	  PasCu = Abs(PasC*previousd2d.X());
+	  PasCv = Abs(PasC*previousd2d.Y());
+
+	  if (PasCu <= tolerance(1) && PasCv <= tolerance(2)) {
+            if (CurrentLine->NbPoints() == 1)  break; // cancel the line
+	    Arrive = Standard_True;
+	    CurrentLine->AddStatusFirstLast(Standard_False,Standard_False,
+                                           Standard_False);
+	    Tgtend = Standard_True;
+	    Rajout = Standard_True;
+	    seqAjout.Append(lines.Length()+1);
+	  }  
+	}
+      }// end of started line 
+      if (Arrive) {
+	CurrentLine->SetTangencyAtBegining(Tgtbeg);
+	CurrentLine->SetTangencyAtEnd(Tgtend);
+        CurrentLine->SetTangency(IsTangentialIntersection);
+	
+	lines.Append(CurrentLine);
+	wd3[I].etat = -wd3[I].etat; //mark point as processed
+      }
+    } //end of processing of tangent point
+  } //end of all tangent points
 }
+
+
+
+
+
+

src/IntWalk/IntWalk_IWalking_5.gxx

@@ -20,6 +20,7 @@ namespace {
                                                    // because the angle is too great in 2d (U4)
 }
 
+
 IntWalk_StatusDeflection IntWalk_IWalking::TestDeflection
   (TheIWFunction& sp,
    const Standard_Boolean Finished,
@@ -27,7 +28,10 @@ IntWalk_StatusDeflection IntWalk_IWalking::TestDeflection
    const IntWalk_StatusDeflection StatusPrecedent,
    Standard_Integer& NbDivision,
    Standard_Real& Step,
-   const Standard_Integer StepSign)
+   const Standard_Integer StepSign,
+   const Standard_Integer CurNbPoints,
+   const Standard_Boolean ArretAjout,
+   const Standard_Boolean OnPrevTangency)
 {
   // Check the step of advancement, AND recalculate this step :
   //
@@ -58,83 +62,90 @@ IntWalk_StatusDeflection IntWalk_IWalking::TestDeflection
   //-- if epsilon is great enough (1e-11). In this case one loops 
   //-- without ever changing the values sent to Rsnld. 
   //---------------------------------------------------------------------------------
-  Standard_Real Paramu = 0.0, Paramv = 0.0;
+  Standard_Real Paramu = 0., Paramv = 0., StepU = 0., StepV = 0.;
-
+  Standard_Real Cosi = 0., Cosi2 = 0., Norme = 0.;
-  if (!reversed) {
-    previousPoint.ParametersOnS2(Paramu, Paramv);
-  }
-  else
-  {
-    previousPoint.ParametersOnS1(Paramu, Paramv);
-  }
-
-  const Standard_Real Du = UV(1) - Paramu;
-  const Standard_Real Dv = UV(2) - Paramv;
-  const Standard_Real Duv = Du * Du + Dv * Dv;
 
   gp_Vec Corde(previousPoint.Value(), sp.Point());
 
-  const Standard_Real Norme = Corde.SquareMagnitude(), 
+  Norme = Corde.SquareMagnitude();
-                      aTol = epsilon*Precision::PConfusion();
+//  if (Norme <= epsilon*epsilon) {
-
+//  if (Norme <= epsilon) { // the square is already taken in the constructor
-  //if ((++NbPointsConfondusConsecutifs < 10) && (Norme <= epsilon)) { // the square is already taken in the constructor
+  Standard_Real TolNorm = gp::Resolution();
-  if ((Norme <= epsilon) && ((Duv <= aTol) || (StatusPrecedent != IntWalk_OK)))
+  if (CurNbPoints == 1 && ArretAjout && OnPrevTangency)
-  { // the square is already taken in the constructor
+    TolNorm = 1.e-10;
+  if (Norme <= TolNorm) {
     Status = IntWalk_PointConfondu;
     if (StatusPrecedent == IntWalk_PasTropGrand) {
       return IntWalk_ArretSurPointPrecedent;
     }
   }
   else {
-    Standard_Real Cosi = Corde * previousd3d;
+    Cosi = Corde * previousd3d;
-    Standard_Real Cosi2 = 0.0;
+    if (Cosi*StepSign < 0.) { // angle 3d > pi/2 !!!!
-
+      Cosi2 = 0.;
-    if (Cosi*StepSign >= 0.) {// angle 3d <= pi/2 !!!!
+    }
-      const Standard_Real aDiv = previousd3d.SquareMagnitude()*Norme;
+    else {
-      if(aDiv == 0)
+      Cosi2 = Cosi * Cosi / previousd3d.SquareMagnitude() / Norme;
-        return Status;
-      Cosi2 = Cosi * Cosi / aDiv;
     }
     if (Cosi2 < CosRef3D) { //angle 3d too great
       Step = Step /2.0;
-      Standard_Real StepU = Abs(Step*previousd2d.X()),
+      StepU = Abs(Step*previousd2d.X());
-                    StepV = Abs(Step*previousd2d.Y());
+      StepV = Abs(Step*previousd2d.Y());
       if (StepU < tolerance(1) && StepV < tolerance(2)) 
-        Status = IntWalk_ArretSurPointPrecedent;
+	Status = IntWalk_ArretSurPointPrecedent;
       else 
-        Status = IntWalk_PasTropGrand;
+	Status = IntWalk_PasTropGrand;
       return Status;
     }
   }
 
-  if (Abs(Du) < tolerance(1) && Abs(Dv) < tolerance(2))
+  if (!reversed) {
+    previousPoint.ParametersOnS2(Paramu, Paramv);
+  }
+  else {
+    previousPoint.ParametersOnS1(Paramu, Paramv);
+  }
+  Standard_Real Du = UV(1) - Paramu;
+  Standard_Real Dv = UV(2) - Paramv;
+  Standard_Real Duv = Du * Du + Dv * Dv;
+  Standard_Real TolDu = tolerance(1), TolDv = tolerance(2);
+  if (CurNbPoints == 1 && ArretAjout && OnPrevTangency)
+    TolDu = TolDv = 1.e-7;
+  if (Abs(Du) < TolDu && Abs(Dv) < TolDv)
     return IntWalk_ArretSurPointPrecedent; //confused point 2d
-
+  Cosi = StepSign * (Du * previousd2d.X() + 
-  Standard_Real Cosi = StepSign * (Du * previousd2d.X() + Dv * previousd2d.Y());
+                     Dv * previousd2d.Y());
-
   if (Cosi < 0 && Status == IntWalk_PointConfondu) 
     return IntWalk_ArretSurPointPrecedent; // leave as step back  
                                            // with confused point
 
-  if (sp.IsTangent()) 
+  if (IsTangentialIntersection)
-    return IntWalk_ArretSurPoint;       
+  {
+    if (sp.IsTangentSmooth()) 
+      return IntWalk_ArretSurPoint;       
+  }
+  else
+  {
+    if (sp.IsTangent()) 
+      return IntWalk_ArretSurPoint;
+  }
 
 //if during routing one has subdivided more than  MaxDivision for each
 //previous step, bug on the square; do nothing (experience U4)
 
-  if ((NbDivision < MaxDivision) && (Status != IntWalk_PointConfondu) && 
+  if (NbDivision < MaxDivision && 
-    (StatusPrecedent!= IntWalk_PointConfondu))
+      Status != IntWalk_PointConfondu && 
-  {
+      StatusPrecedent!= IntWalk_PointConfondu ) {
-    Standard_Real Cosi2 = Cosi * Cosi / Duv;
+    Cosi2 = Cosi * Cosi / Duv;
     if (Cosi2 < CosRef2D || Cosi < 0  ) {
       Step = Step / 2.0;
-      Standard_Real StepU = Abs(Step*previousd2d.X()),
+      StepU = Abs(Step*previousd2d.X());
-                    StepV = Abs(Step*previousd2d.Y());
+      StepV = Abs(Step*previousd2d.Y());
 
       if (StepU < tolerance(1) && StepV < tolerance(2))
-        Status = IntWalk_ArretSurPointPrecedent;
+	Status = IntWalk_ArretSurPointPrecedent;
       else 
-        Status = IntWalk_PasTropGrand;
+	Status = IntWalk_PasTropGrand;
       NbDivision = NbDivision + 1;
       return Status;
     }
@@ -143,149 +154,136 @@ IntWalk_StatusDeflection IntWalk_IWalking::TestDeflection
     Cosi2 = Cosi * Cosi / sp.Direction3d().SquareMagnitude() / Norme;
     if (Cosi2 < CosRef3D ){ //angle 3d too great
       Step = Step / 2.;
-      Standard_Real StepU = Abs(Step*previousd2d.X()),
+      StepU = Abs(Step*previousd2d.X());
-                    StepV = Abs(Step*previousd2d.Y());
+      StepV = Abs(Step*previousd2d.Y());
       if (StepU < tolerance(1) && StepV < tolerance(2))
-        Status = IntWalk_ArretSurPoint;
+	Status = IntWalk_ArretSurPoint;
       else 
-        Status = IntWalk_PasTropGrand;
+	Status = IntWalk_PasTropGrand;
       return Status;
     }
     Cosi = Du * sp.Direction2d().X() + 
-      Dv * sp.Direction2d().Y();
+           Dv * sp.Direction2d().Y();
     Cosi2 = Cosi * Cosi / Duv;
     if (Cosi2 < CosRef2D || 
-      sp.Direction2d() * previousd2d < 0) {
+	sp.Direction2d() * previousd2d < 0) {
-        //angle 2d too great or change the side       
+      //angle 2d too great or change the side       
-        Step  = Step / 2.;
+      Step  = Step / 2.;
-        Standard_Real StepU = Abs(Step*previousd2d.X()),
+      StepU = Abs(Step*previousd2d.X());
-                      StepV = Abs(Step*previousd2d.Y());
+      StepV = Abs(Step*previousd2d.Y());
-        if (StepU < tolerance(1) && StepV < tolerance(2))
+      if (StepU < tolerance(1) && StepV < tolerance(2))
-          Status = IntWalk_ArretSurPointPrecedent;
+	Status = IntWalk_ArretSurPointPrecedent;
-        else 
+      else 
-          Status = IntWalk_PasTropGrand;
+	Status = IntWalk_PasTropGrand;
-        return Status;
+      return Status;
     }
   }
 
   if (!Finished) {
-    if (Status == IntWalk_PointConfondu)
+
-    {
+    if (Status == IntWalk_PointConfondu) {
-      Standard_Real StepU = Min(Abs(1.5 * Du),pas*(UM-Um)),
+      StepU = Min(Abs(1.5 * Du),pas*(UM-Um));
-                    StepV = Min(Abs(1.5 * Dv),pas*(VM-Vm));
+      StepV = Min(Abs(1.5 * Dv),pas*(VM-Vm));
 
       Standard_Real d2dx = Abs(previousd2d.X()); 
       Standard_Real d2dy = Abs(previousd2d.Y()); 
 
-      if (d2dx < tolerance(1))
+      if (d2dx < tolerance(1)) {
-      {
+	Step = StepV/d2dy;
-        Step = StepV/d2dy;
       }
-      else if (d2dy < tolerance(2))
+      else if (d2dy < tolerance(2)) {
-      {
+	Step = StepU/d2dx;
-        Step = StepU/d2dx;
       }
-      else
+      else {
-      {
+	Step = Min(StepU/d2dx,StepV/d2dy);
-        Step = Min(StepU/d2dx,StepV/d2dy);
       }
     }
-    else
+    else {
-    {
+//   estimate the current vector.
-      //   estimate the current vector.
+//   if vector/2<=current vector<= vector it is considered that the criterion
-      //   if vector/2<=current vector<= vector it is considered that the criterion
+//   is observed.
-      //   is observed.
+//   otherwise adjust the step depending on the previous step 
-      //   otherwise adjust the step depending on the previous step 
 
-      /*
+/*
-        Standard_Real Dist = Sqrt(Norme)/3.;
+      Standard_Real Dist = Sqrt(Norme)/3.;
-        TColgp_Array1OfPnt Poles(1,4);
+      TColgp_Array1OfPnt Poles(1,4);
-        gp_Pnt POnCurv,Milieu;
+      gp_Pnt POnCurv,Milieu;
-        Poles(1) = previousPoint.Value();
+      Poles(1) = previousPoint.Value();
-        Poles(4) = sp.Point();
+      Poles(4) = sp.Point();
-        Poles(2) = Poles(1).XYZ() + 
+      Poles(2) = Poles(1).XYZ() + 
-      StepSign * Dist* previousd3d.Normalized().XYZ();
+	StepSign * Dist* previousd3d.Normalized().XYZ();
-        Poles(3) = Poles(4).XYZ() - 
+      Poles(3) = Poles(4).XYZ() - 
-      StepSign * Dist*sp.Direction3d().Normalized().XYZ();
+	StepSign * Dist*sp.Direction3d().Normalized().XYZ();
-        BzCLib::PntPole(0.5,Poles,POnCurv);
+      BzCLib::PntPole(0.5,Poles,POnCurv);
-        Milieu = (Poles(1).XYZ() + Poles(4).XYZ())*0.5;
+      Milieu = (Poles(1).XYZ() + Poles(4).XYZ())*0.5;
-      //      FlecheCourante = Milieu.Distance(POnCurv);
+//      FlecheCourante = Milieu.Distance(POnCurv);
-        Standard_Real FlecheCourante = Milieu.SquareDistance(POnCurv);
+      Standard_Real FlecheCourante = Milieu.SquareDistance(POnCurv);
-      */
+*/
 
-        // Direct calculation : 
+      // Direct calculation : 
-        // POnCurv=(((p1+p2)/2.+(p2+p3)/2.)/2. + ((p2+p3)/2.+(p3+P4)/2.)/2.)/2.
+      // POnCurv=(((p1+p2)/2.+(p2+p3)/2.)/2. + ((p2+p3)/2.+(p3+P4)/2.)/2.)/2.
-        // either POnCurv = p1/8. + 3.p2/8. + 3.p3/8. + p4/8.
+      // either POnCurv = p1/8. + 3.p2/8. + 3.p3/8. + p4/8.
-        // Or p2 = p1 + lambda*d1 et p3 = p4 - lambda*d4
+      // Or p2 = p1 + lambda*d1 et p3 = p4 - lambda*d4
-        // So POnCurv = (p1 + p4)/2. + 3.*(lambda d1 - lambda d4)/8.
+      // So POnCurv = (p1 + p4)/2. + 3.*(lambda d1 - lambda d4)/8.
-        // Calculate the deviation with (p1+p4)/2. . So it is just necessary to calculate
+      // Calculate the deviation with (p1+p4)/2. . So it is just necessary to calculate
-        // the norm (square) of 3.*lambda (d1 - d4)/8.
+      // the norm (square) of 3.*lambda (d1 - d4)/8.
-        // either the norm of :
+      // either the norm of :
-        //    3.*(Sqrt(Norme)/3.)*StepSign*(d1-d4)/8.
+      //    3.*(Sqrt(Norme)/3.)*StepSign*(d1-d4)/8.
-        // which produces, takin the square :
+      // which produces, takin the square :
-        //         Norme * (d1-d4).SquareMagnitude()/64.
+      //         Norme * (d1-d4).SquareMagnitude()/64.
 
       Standard_Real FlecheCourante = 
 	(previousd3d.Normalized().XYZ()-sp.Direction3d().Normalized().XYZ()).SquareModulus()*Norme/64.;
 
   
 //      if (FlecheCourante <= 0.5*fleche) {
-      if (FlecheCourante <= 0.25*fleche*fleche)
+      if (FlecheCourante <= 0.25*fleche*fleche) {
-      {
-        Standard_Real d2dx = Abs(sp.Direction2d().X()); 
-        Standard_Real d2dy = Abs(sp.Direction2d().Y()); 
-        
-        Standard_Real StepU = Min(Abs(1.5*Du),pas*(UM-Um)),
-                      StepV = Min(Abs(1.5*Dv),pas*(VM-Vm));
 
-        if (d2dx < tolerance(1))
+	Standard_Real d2dx = Abs(sp.Direction2d().X()); 
-        {
+	Standard_Real d2dy = Abs(sp.Direction2d().Y()); 
-          Step = StepV/d2dy;
+
-        }
+	StepU = Min(Abs(1.5*Du),pas*(UM-Um));
-        else if (d2dy < tolerance(2))
+	StepV = Min(Abs(1.5*Dv),pas*(VM-Vm));
-        {
+
-          Step = StepU/d2dx;
+	if (d2dx < tolerance(1)) {
-        }
+	  Step = StepV/d2dy;
-        else
+	}
-        {
+	else if (d2dy < tolerance(2)) {
-          Step = Min(StepU/d2dx,StepV/d2dy);
+	  Step = StepU/d2dx;
-        }	
+	}
+	else {
+	  Step = Min(StepU/d2dx,StepV/d2dy);
+	}
+	
       }
-      else
+      else {
-      {
+//	if (FlecheCourante > fleche) {  // step too great
-        //if (FlecheCourante > fleche) {  // step too great
+	if (FlecheCourante > fleche*fleche) {  // step too great
-        if (FlecheCourante > fleche*fleche)
+	  Step = Step /2.;
-        {  // step too great
+          StepU = Abs(Step*previousd2d.X());
-          Step = Step /2.;
+          StepV = Abs(Step*previousd2d.Y());
-          Standard_Real StepU = Abs(Step*previousd2d.X()),
-                        StepV = Abs(Step*previousd2d.Y());
-          
           if (StepU < tolerance(1) && StepV < tolerance(2)) 
             Status = IntWalk_ArretSurPointPrecedent;
           else 
             Status = IntWalk_PasTropGrand;
-        }
+	}
-        else
+	else {
-        {
+	  Standard_Real d2dx = Abs(sp.Direction2d().X()); 
-          Standard_Real d2dx = Abs(sp.Direction2d().X()); 
+	  Standard_Real d2dy = Abs(sp.Direction2d().Y()); 
-          Standard_Real d2dy = Abs(sp.Direction2d().Y()); 
-          
-          Standard_Real StepU = Min(Abs(1.5*Du),pas*(UM-Um)),
-                        StepV = Min(Abs(1.5*Dv),pas*(VM-Vm));
 
-          if (d2dx < tolerance(1))
+	  StepU = Min(Abs(1.5*Du),pas*(UM-Um));
-          {
+	  StepV = Min(Abs(1.5*Dv),pas*(VM-Vm));
-            Step = Min(Step,StepV/d2dy);
+
-          }
+	  if (d2dx < tolerance(1)) {
-          else if (d2dy < tolerance(2))
+	    Step = Min(Step,StepV/d2dy);
-          {
+	  }
-            Step = Min(Step,StepU/d2dx);
+	  else if (d2dy < tolerance(2)) {
-          }
+	    Step = Min(Step,StepU/d2dx);
-          else
+	  }
-          {
+	  else {
-            Step = Min(Step,Min(StepU/d2dx,StepV/d2dy));
+	    Step = Min(Step,Min(StepU/d2dx,StepV/d2dy));
-          }
+	  }
-        }
+	}
       }
     }
   }

src/IntWalk/IntWalk_IWalking_6.gxx

@@ -12,12 +12,112 @@
 // Alternatively, this file may be used under the terms of Open CASCADE
 // commercial license or contractual agreement.
 
+#include <Geom2d_Line.hxx>
+#include <GCE2d_MakeLine.hxx>
+
+
 #ifndef OCCT_DEBUG
 #define No_Standard_RangeError
 #define No_Standard_OutOfRange
 #endif
 
 
+const Standard_Real TolDiscr = 1.e-8;
+const Standard_Real TolB = 1.e-5;
+
+static Standard_Integer RealSign(const Standard_Real aValue)
+{
+  if (aValue >= 0.)
+    return 1;
+  else return -1;
+}
+
+static Standard_Boolean ComputeB11orB22(TheIWFunction& sp,
+                                        Standard_Boolean& IsDiscriminantNull,
+                                        Standard_Integer& SignOfBcoeff,
+                                        gp_Vec& Pu,
+                                        gp_Vec& Pv,
+                                        Standard_Real& B11,
+                                        Standard_Real& B22,
+                                        Standard_Real& B12)
+{
+  gp_Vec Puu, Puv, Pvv;
+  gp_Vec Iu, Iv, Iuu, Iuv, Ivv;
+  gp_Vec NormalP, NormalI;
+
+  if (!sp.DerivativesAndNormalOnPSurf(Pu, Pv, NormalP, Puu, Pvv, Puv))
+    return Standard_False;
+  if (!sp.DerivativesAndNormalOnISurf(Iu, Iv, NormalI, Iuu, Ivv, Iuv))
+    return Standard_False;
+
+  Standard_Real Lp, Mp, Np, Li, Mi, Ni; //second fundamental form coefficients
+  Lp = Puu * NormalP;
+  Mp = Puv * NormalP;
+  Np = Pvv * NormalP;
+  Li = Iuu * NormalI;
+  Mi = Iuv * NormalI;
+  Ni = Ivv * NormalI;
+
+  gp_Vec Normal;
+  if (NormalP * NormalI < 0.)
+    NormalI.Reverse();
+  Normal.SetXYZ(0.5*(NormalP.XYZ() + NormalI.XYZ()));
+  Normal.Normalize();
+
+  Standard_Real A11, A12, A21, A22;
+  Standard_Real NormIdotNorm = (Iu ^ Iv) * Normal;
+  A11 = ((Pu ^ Iv) * Normal) / NormIdotNorm;
+  A12 = ((Pv ^ Iv) * Normal) / NormIdotNorm;
+  A21 = ((Iu ^ Pu) * Normal) / NormIdotNorm;
+  A22 = ((Iu ^ Pv) * Normal) / NormIdotNorm;
+
+  //Standard_Real B11, B12, B22;
+  B11 = A11*A11*Li + 2*A11*A21*Mi + A21*A21*Ni - Lp;
+  B12 = A11*A12*Li + (A11*A22 + A21*A12)*Mi + A21*A22*Ni - Mp;
+  B22 = A12*A12*Li + 2*A12*A22*Mi + A22*A22*Ni - Np;
+
+  Standard_Real Discriminant = B12*B12 - B11*B22;
+
+  if (Abs(Discriminant) < TolDiscr)
+  {
+    IsDiscriminantNull = Standard_True;
+    if (Abs(B11) <= TolB && Abs(B22) <= TolB)
+    {
+#ifdef DRAW
+      cout<<"Possible branching"<<endl<<endl;
+#endif
+      SignOfBcoeff = 0;
+    }
+    else
+    {
+      gp_Vec NewPreviousD3d;
+      gp_Dir2d NewPreviousD2d;
+      if (Abs(B11) > TolB)
+      {
+#ifdef DRAW
+        cout<<"B11 = "<<B11<<endl;
+#endif
+        SignOfBcoeff = RealSign(B11);
+      }
+      else
+      {
+#ifdef DRAW
+        cout<<"B22 = "<<B22<<endl;
+#endif
+        SignOfBcoeff = RealSign(B22);
+      }
+    }
+  }
+  else
+  {
+    IsDiscriminantNull = Standard_False;
+#ifdef DRAW
+    cout<<"Discriminant = "<<Discriminant<<endl;
+#endif
+  }
+  return Standard_True;
+}
+
 void IntWalk_IWalking::AddPointInCurrentLine
          (const Standard_Integer N,
 	  const ThePointOfPath& PathPnt,
@@ -39,7 +139,6 @@ void IntWalk_IWalking::MakeWalkingPoint
 	  IntSurf_PntOn2S& Psol )
 
 {
-
 // Case == 1      : make a WalkinkPoint.
 // Case == 2      : make a WalkinkPoint.
 //                  The computation of the tangency on is done  
@@ -60,15 +159,117 @@ void IntWalk_IWalking::MakeWalkingPoint
     UV(1) = U;
     UV(2) = V;
     sp.Values(UV, FF, DD);
-    MakeWalkingPoint(Case - 10, U, V, sp, Psol); 
+    MakeWalkingPoint(Case - 10, U, V, sp, Psol);
   }
   else {
     Standard_ConstructionError::Raise();
   }
-
-
 }
 
+Standard_Boolean
+ IntWalk_IWalking::ComputeDirOfTangentialIntersection(TheIWFunction& sp,
+                                                      Standard_Integer& StepSignTangent,
+                                                      Standard_Boolean& IsDiscriminantNull,
+                                                      Standard_Integer& SignOfBcoeff,
+                                                      gp_Vec&   newd3d,
+                                                      gp_Dir2d& newd2d)
+{
+  gp_Vec Pu, Pv;
+  Standard_Real B11, B22, B12;
+  if (!ComputeB11orB22(sp, IsDiscriminantNull, SignOfBcoeff,
+                       Pu, Pv, B11, B22, B12))
+    return Standard_False;
+  
+  if (IsDiscriminantNull)
+  {
+    if (SignOfBcoeff != 0)
+    {
+      gp_Vec NewPreviousD3d;
+      gp_Dir2d NewPreviousD2d;
+      if (Abs(B11) > TolB)
+      {
+#ifdef DRAW
+        cout<<"B11 = "<<B11<<endl;
+#endif
+        Standard_Real CoefPu = -B12/B11;
+        NewPreviousD3d = CoefPu*Pu + Pv;
+        NewPreviousD3d.Normalize();
+        NewPreviousD2d = gp_Dir2d(CoefPu, 1.);
+      }
+      else
+      {
+#ifdef DRAW
+        cout<<"B22 = "<<B22<<endl;
+#endif
+        Standard_Real CoefPv = -B12/B22;
+        NewPreviousD3d = Pu + CoefPv*Pv;
+        NewPreviousD3d.Normalize();
+        NewPreviousD2d = gp_Dir2d(1., CoefPv);
+      }
+      
+      if (!IsTangentialIntersection)
+      {
+        IsTangentialIntersection = Standard_True;
+        if (NewPreviousD3d * previousd3d < 0)
+          StepSignTangent = -1;
+        else
+          StepSignTangent = 1;
+      }
+      newd3d = StepSignTangent * NewPreviousD3d;
+      newd2d = StepSignTangent * NewPreviousD2d;
+    }
+  }
+  return Standard_True;
+}
+
+void IntWalk_IWalking::FindExactCuspPoint(TheIWFunction& sp,
+                                          IntSurf_PntOn2S& Psol,
+                                          const Standard_Integer SignFirst,
+                                          const Standard_Integer SignLast)
+{
+  IntSurf_PntOn2S EndPoint = Psol;
+
+  IntSurf_PntOn2S Pfirst, Plast;
+
+  Standard_Real newU, newV;
+  IntSurf_PntOn2S newPoint;
+
+  Pfirst = previousPoint;
+  Plast =  EndPoint;
+
+  Standard_Integer newSign;
+
+  for (;;)
+  {
+    Standard_Real Ufirst, Vfirst, Ulast, Vlast;
+    Pfirst.ParametersOnSurface(reversed, Ufirst, Vfirst);
+    Plast.ParametersOnSurface(reversed, Ulast, Vlast);
+    
+    newU = 0.5*(Ufirst + Ulast);
+    newV = 0.5*(Vfirst + Vlast);
+
+    if (Abs(newU - Ufirst) < tolerance(1) &&
+        Abs(newV - Vfirst) < tolerance(2))
+      break;
+    
+    MakeWalkingPoint(11, newU, newV, sp, newPoint);
+    Standard_Boolean IsDiscriminantNull;
+    gp_Vec Pu, Pv;
+    Standard_Real B11, B22, B12;
+    ComputeB11orB22(sp, IsDiscriminantNull, newSign,
+                    Pu, Pv, B11, B22, B12);
+    
+    if (newSign == 0)
+      break;
+
+    if (newSign == SignFirst)
+      Pfirst = newPoint;
+    else //newSign == SignLast
+      Plast = newPoint;
+  }
+
+  Psol = newPoint;
+}
 
 
 void IntWalk_IWalking::OpenLine(const Standard_Integer N,
@@ -94,8 +295,26 @@ void IntWalk_IWalking::OpenLine(const Standard_Integer N,
     previousPoint.ParametersOnS1(UV(1),UV(2));
   }
   sp.Values(UV, FF, DD);
-  previousd3d = sp.Direction3d();
+  if (sp.IsTangent())
-  previousd2d = sp.Direction2d();
+  {
+    IsTangentialIntersection = Standard_True;
+    Standard_Integer theSign = 1;
+    Standard_Integer theSignOfBcoeff = 0;
+    gp_Vec   newd3d;
+    gp_Dir2d newd2d;
+    Standard_Boolean IsDiscriminantNull;
+    ComputeDirOfTangentialIntersection(sp, theSign,
+                                       IsDiscriminantNull, theSignOfBcoeff,
+                                       newd3d, newd2d);
+    previousd3d = newd3d;
+    previousd2d = newd2d;
+  }
+  else
+  {
+    IsTangentialIntersection = Standard_False;
+    previousd3d = sp.Direction3d();
+    previousd2d = sp.Direction2d();
+  }
 
   if (N>0) { //departure point given at input
     PathPnt = Pnts1.Value(N);