0030675: Visualization - remove redundant proxy classes in hierarchy of PrsMgr_Presen...
[occt.git] / src / IntTools / IntTools_FaceFace.cxx
index 40cebc2..759e5dd 100644 (file)
 // Alternatively, this file may be used under the terms of Open CASCADE
 // commercial license or contractual agreement.
 
-#include <IntTools_FaceFace.ixx>
+#include <IntTools_FaceFace.hxx>
 
-#include <Precision.hxx>
-
-#include <TColStd_HArray1OfReal.hxx>
-#include <TColStd_Array1OfReal.hxx>
-#include <TColStd_Array1OfInteger.hxx>
-#include <TColStd_SequenceOfReal.hxx>
-#include <TColStd_ListOfInteger.hxx>
-#include <TColStd_ListIteratorOfListOfInteger.hxx>
-#include <TColStd_Array1OfListOfInteger.hxx>
-
-#include <gp_Lin2d.hxx>
-#include <gp_Ax22d.hxx>
-#include <gp_Circ2d.hxx>
-#include <gp_Torus.hxx>
-#include <gp_Cylinder.hxx>
-
-#include <Bnd_Box.hxx>
-
-#include <TColgp_HArray1OfPnt2d.hxx>
-#include <TColgp_SequenceOfPnt2d.hxx>
-#include <TColgp_Array1OfPnt.hxx>
-#include <TColgp_Array1OfPnt2d.hxx>
-
-#include <IntAna_QuadQuadGeo.hxx>
-
-#include <IntSurf_PntOn2S.hxx>
-#include <IntSurf_LineOn2S.hxx>
-#include <IntSurf_PntOn2S.hxx>
-#include <IntSurf_ListOfPntOn2S.hxx>
-#include <IntRes2d_Domain.hxx>
-#include <ProjLib_Plane.hxx>
-
-#include <IntPatch_GLine.hxx>
-#include <IntPatch_RLine.hxx>
-#include <IntPatch_WLine.hxx>
-#include <IntPatch_ALine.hxx>
-#include <IntPatch_ALineToWLine.hxx>
-
-#include <ElSLib.hxx>
+#include <BRepAdaptor_Surface.hxx>
+#include <BRepTools.hxx>
+#include <BRep_Tool.hxx>
 #include <ElCLib.hxx>
-
-#include <Extrema_ExtCC.hxx>
-#include <Extrema_POnCurv.hxx>
-#include <BndLib_AddSurface.hxx>
-
-#include <Adaptor3d_SurfacePtr.hxx>
-#include <Adaptor2d_HLine2d.hxx>
-
-#include <GeomAbs_SurfaceType.hxx>
-#include <GeomAbs_CurveType.hxx>
-
-#include <Geom_Surface.hxx>
-#include <Geom_Line.hxx>
+#include <ElSLib.hxx>
+#include <Geom2dAdaptor_Curve.hxx>
+#include <Geom2dInt_GInter.hxx>
+#include <Geom2d_Line.hxx>
+#include <Geom2d_TrimmedCurve.hxx>
+#include <GeomAPI_ProjectPointOnSurf.hxx>
+#include <GeomAdaptor_HSurface.hxx>
+#include <GeomInt_IntSS.hxx>
+#include <GeomInt_WLApprox.hxx>
+#include <GeomLib_Check2dBSplineCurve.hxx>
+#include <GeomLib_CheckBSplineCurve.hxx>
+#include <Geom_BSplineCurve.hxx>
 #include <Geom_Circle.hxx>
 #include <Geom_Ellipse.hxx>
-#include <Geom_Parabola.hxx>
 #include <Geom_Hyperbola.hxx>
-#include <Geom_TrimmedCurve.hxx>
-#include <Geom_BSplineCurve.hxx>
-#include <Geom_RectangularTrimmedSurface.hxx>
+#include <Geom_Line.hxx>
 #include <Geom_OffsetSurface.hxx>
-#include <Geom_Curve.hxx>
-#include <Geom_Conic.hxx>
-
-#include <Geom2d_TrimmedCurve.hxx>
-#include <Geom2d_BSplineCurve.hxx>
-#include <Geom2d_Line.hxx>
-#include <Geom2d_Curve.hxx>
-#include <Geom2d_Circle.hxx>
-
-#include <Geom2dAPI_InterCurveCurve.hxx>
-#include <Geom2dInt_GInter.hxx>
-#include <Geom2dAdaptor.hxx>
-#include <GeomAdaptor_Curve.hxx>
-#include <GeomAdaptor_HSurface.hxx>
-#include <GeomAdaptor_Surface.hxx>
-#include <GeomLib_CheckBSplineCurve.hxx>
-#include <GeomLib_Check2dBSplineCurve.hxx>
-
-#include <GeomInt_WLApprox.hxx>
-#include <GeomProjLib.hxx>
-#include <GeomAPI_ProjectPointOnSurf.hxx>
-#include <Geom2dAdaptor_Curve.hxx>
-#include <TopoDS.hxx>
-#include <TopoDS_Edge.hxx>
-#include <TopExp_Explorer.hxx>
-
-#include <BRep_Tool.hxx>
-#include <BRepTools.hxx>
-#include <BRepAdaptor_Surface.hxx>
-
-#include <IntTools_Curve.hxx>
-#include <IntTools_Tools.hxx>
+#include <Geom_Parabola.hxx>
+#include <Geom_RectangularTrimmedSurface.hxx>
+#include <Geom_TrimmedCurve.hxx>
+#include <IntAna_QuadQuadGeo.hxx>
+#include <IntPatch_GLine.hxx>
+#include <IntPatch_RLine.hxx>
+#include <IntPatch_WLine.hxx>
+#include <IntRes2d_Domain.hxx>
+#include <IntSurf_Quadric.hxx>
+#include <IntTools_Context.hxx>
 #include <IntTools_Tools.hxx>
 #include <IntTools_TopolTool.hxx>
-#include <IntTools_PntOnFace.hxx>
-#include <IntTools_PntOn2Faces.hxx>
-#include <IntTools_Context.hxx>
-#include <IntSurf_ListIteratorOfListOfPntOn2S.hxx>
-#include <GeomInt.hxx>
-
-#include <Approx_CurveOnSurface.hxx>
-#include <GeomAdaptor.hxx>
-#include <GeomInt_IntSS.hxx>
-
-static
-  void RefineVector(gp_Vec2d& aV2D);
-#ifdef OCCT_DEBUG_DUMPWLINE
-static
-  void DumpWLine(const Handle(IntPatch_WLine)& aWLine);
-#endif
-//
-static
-  void TolR3d(const TopoDS_Face& ,
-              const TopoDS_Face& ,
-              Standard_Real& );
-static 
-  Handle(Geom_Curve) MakeBSpline (const Handle(IntPatch_WLine)&,
-                                  const Standard_Integer,
-                                  const Standard_Integer);
+#include <IntTools_WLineTool.hxx>
+#include <ProjLib_Plane.hxx>
+#include <TopExp_Explorer.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Edge.hxx>
+#include <gp_Elips.hxx>
 
 static 
   void Parameters(const Handle(GeomAdaptor_HSurface)&,
@@ -153,35 +72,6 @@ static
                                 Standard_Real&      thevmin, 
                                 Standard_Real&      thevmax);
 
-static
-  Standard_Boolean NotUseSurfacesForApprox
-          (const TopoDS_Face& aF1,
-           const TopoDS_Face& aF2,
-           const Handle(IntPatch_WLine)& WL,
-           const Standard_Integer ifprm,
-           const Standard_Integer ilprm);
-
-static 
-  Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine);
-
-static 
-  Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
-                                            const Standard_Integer                         ideb,
-                                            const Standard_Integer                         ifin,
-                                            const Standard_Boolean                         onFirst);
-
-static 
-  Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
-                                        const Handle(GeomAdaptor_HSurface)&            theSurface1, 
-                                        const Handle(GeomAdaptor_HSurface)&            theSurface2,
-                                        const TopoDS_Face&                             theFace1,
-                                        const TopoDS_Face&                             theFace2,
-                                        const GeomInt_LineConstructor&                 theLConstructor,
-                                        const Standard_Boolean                         theAvoidLConstructor,
-                                        IntPatch_SequenceOfLine&                       theNewLines,
-                                        Standard_Real&                                 theReachedTol3d,
-                                        const Handle(IntTools_Context)& );
-
 static 
   Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter, 
                                           const Handle(Geom_Curve)& theCurve, 
@@ -189,71 +79,26 @@ static
                                           const TopoDS_Face&        theFace2,
                                           const Standard_Real       theOtherParameter,
                                           const Standard_Boolean    bIncreasePar,
+                                          const Standard_Real       theTol,
                                           Standard_Real&            theNewParameter,
                                           const Handle(IntTools_Context)& );
 
 static 
-  Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve);
-
-static 
-  Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
-                                     const Standard_Real theFirstBoundary,
-                                     const Standard_Real theSecondBoundary,
-                                     const Standard_Real theResolution,
-                                     Standard_Boolean&   IsOnFirstBoundary);
-static
-  Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
-                             const gp_Pnt2d&     theLastPoint,
-                             const Standard_Real theUmin, 
-                             const Standard_Real theUmax,
-                             const Standard_Real theVmin,
-                             const Standard_Real theVmax,
-                             gp_Pnt2d&           theNewPoint);
-
-
-static 
-  Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
-                                       const Handle(GeomAdaptor_HSurface)&  theSurface2,
-                                       const TopoDS_Face&                   theFace1,
-                                       const TopoDS_Face&                   theFace2,
-                                       Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
-                                       Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
-                                       Handle(TColStd_HArray1OfReal)&       theResultRadius,
-                                       const Handle(IntTools_Context)& );
-
-static
-  Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
-                             const gp_Pnt2d&     theLastPoint,
-                             const Standard_Real theUmin, 
-                             const Standard_Real theUmax,
-                             const Standard_Real theVmin,
-                             const Standard_Real theVmax,
-                             const gp_Pnt2d&     theTanZoneCenter,
-                             const Standard_Real theZoneRadius,
-                             Handle(GeomAdaptor_HSurface) theGASurface,
-                             gp_Pnt2d&           theNewPoint);
+  Standard_Boolean IsCurveValid(const Handle(Geom2d_Curve)& thePCurve);
 
-static
-  Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
-                                   const gp_Pnt2d&     theTanZoneCenter,
-                                   const Standard_Real theZoneRadius,
-                                   Handle(GeomAdaptor_HSurface) theGASurface);
-
-static
-  gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
-                             const gp_Pnt2d&     theOriginalPoint,
-                             Handle(GeomAdaptor_HSurface) theGASurface);
 static
   Standard_Boolean  ApproxWithPCurves(const gp_Cylinder& theCyl, 
                                       const gp_Sphere& theSph);
 
-static void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1, 
-                           const Handle(GeomAdaptor_HSurface)& theS2, 
-                           const Standard_Real TolAng, 
-                           const Standard_Real TolTang, 
+static void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1,
+                           const Handle(GeomAdaptor_HSurface)& theS2,
+                           const Standard_Real TolF1,
+                           const Standard_Real TolF2,
+                           const Standard_Real TolAng,
+                           const Standard_Real TolTang,
                            const Standard_Boolean theApprox1,
                            const Standard_Boolean theApprox2,
-                           IntTools_SequenceOfCurves& theSeqOfCurve, 
+                           IntTools_SequenceOfCurves& theSeqOfCurve,
                            Standard_Boolean& theTangentFaces);
 
 static Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS, 
@@ -283,7 +128,27 @@ static
 //
 static
   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC, 
-                               const TopoDS_Face& aFace);
+                               const TopoDS_Face& aFace,
+                               const Handle(IntTools_Context)& theCtx);
+
+static
+  Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
+                            const Standard_Real aT,
+                            GeomAPI_ProjectPointOnSurf& theProjPS);
+
+static
+  Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
+                                const Standard_Real theFirst,
+                                const Standard_Real theLast,
+                                GeomAPI_ProjectPointOnSurf& theProjPS,
+                                const Standard_Real theEps);
+
+static
+  Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
+                                const Standard_Real theFirst,
+                                const Standard_Real theLast,
+                                const TopoDS_Face& theFace,
+                                const Handle(IntTools_Context)& theContext);
 
 static
   void CorrectPlaneBoundaries(Standard_Real& aUmin,
@@ -302,10 +167,11 @@ IntTools_FaceFace::IntTools_FaceFace()
   //
   myHS1 = new GeomAdaptor_HSurface ();
   myHS2 = new GeomAdaptor_HSurface ();
-  myTolReached2d=0.; 
-  myTolReached3d=0.;
+  myTolF1 = 0.;
+  myTolF2 = 0.;
+  myTol = 0.;
+  myFuzzyValue = Precision::Confusion();
   SetParameters(Standard_True, Standard_True, Standard_True, 1.e-07);
-  
 }
 //=======================================================================
 //function : SetContext
@@ -364,14 +230,6 @@ Standard_Boolean IntTools_FaceFace::IsDone() const
   return myIsDone;
 }
 //=======================================================================
-//function : TolReached3d
-//purpose  : 
-//=======================================================================
-Standard_Real IntTools_FaceFace::TolReached3d() const
-{
-  return myTolReached3d;
-}
-//=======================================================================
 //function : Lines
 //purpose  : return lines of intersection
 //=======================================================================
@@ -382,14 +240,6 @@ const IntTools_SequenceOfCurves& IntTools_FaceFace::Lines() const
      "IntTools_FaceFace::Lines() => myIntersector NOT DONE");
   return mySeqOfCurve;
 }
-//=======================================================================
-//function : TolReached2d
-//purpose  : 
-//=======================================================================
-Standard_Real IntTools_FaceFace::TolReached2d() const
-{
-  return myTolReached2d;
-}
 // =======================================================================
 // function: SetParameters
 //
@@ -404,6 +254,23 @@ void IntTools_FaceFace::SetParameters(const Standard_Boolean ToApproxC3d,
   myApprox2 = ToApproxC2dOnS2;
   myTolApprox = ApproximationTolerance;
 }
+//=======================================================================
+//function : SetFuzzyValue
+//purpose  : 
+//=======================================================================
+void IntTools_FaceFace::SetFuzzyValue(const Standard_Real theFuzz)
+{
+  myFuzzyValue = Max(theFuzz, Precision::Confusion());
+}
+//=======================================================================
+//function : FuzzyValue
+//purpose  : 
+//=======================================================================
+Standard_Real IntTools_FaceFace::FuzzyValue() const
+{
+  return myFuzzyValue;
+}
+
 //=======================================================================
 //function : SetList
 //purpose  : 
@@ -414,28 +281,25 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
 }
 
 
-static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
-                                        const TopoDS_Face& theF2)
+static Standard_Boolean isTreatAnalityc(const BRepAdaptor_Surface& theBAS1,
+                                        const BRepAdaptor_Surface& theBAS2,
+                                        const Standard_Real theTol)
 {
   const Standard_Real Tolang = 1.e-8;
-  const Standard_Real aTolF1=BRep_Tool::Tolerance(theF1);
-  const Standard_Real aTolF2=BRep_Tool::Tolerance(theF2);
-  const Standard_Real aTolSum = aTolF1 + aTolF2;
   Standard_Real aHigh = 0.0;
 
-  const BRepAdaptor_Surface aBAS1(theF1), aBAS2(theF2);
-  const GeomAbs_SurfaceType aType1=aBAS1.GetType();
-  const GeomAbs_SurfaceType aType2=aBAS2.GetType();
+  const GeomAbs_SurfaceType aType1=theBAS1.GetType();
+  const GeomAbs_SurfaceType aType2=theBAS2.GetType();
   
   gp_Pln aS1;
   gp_Cylinder aS2;
   if(aType1 == GeomAbs_Plane)
   {
-    aS1=aBAS1.Plane();
+    aS1=theBAS1.Plane();
   }
   else if(aType2 == GeomAbs_Plane)
   {
-    aS1=aBAS2.Plane();
+    aS1=theBAS2.Plane();
   }
   else
   {
@@ -444,9 +308,9 @@ static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
 
   if(aType1 == GeomAbs_Cylinder)
   {
-    aS2=aBAS1.Cylinder();
-    const Standard_Real VMin = aBAS1.FirstVParameter();
-    const Standard_Real VMax = aBAS1.LastVParameter();
+    aS2=theBAS1.Cylinder();
+    const Standard_Real VMin = theBAS1.FirstVParameter();
+    const Standard_Real VMax = theBAS1.LastVParameter();
 
     if( Precision::IsNegativeInfinite(VMin) ||
         Precision::IsPositiveInfinite(VMax))
@@ -456,10 +320,10 @@ static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
   }
   else if(aType2 == GeomAbs_Cylinder)
   {
-    aS2=aBAS2.Cylinder();
+    aS2=theBAS2.Cylinder();
 
-    const Standard_Real VMin = aBAS2.FirstVParameter();
-    const Standard_Real VMax = aBAS2.LastVParameter();
+    const Standard_Real VMin = theBAS2.FirstVParameter();
+    const Standard_Real VMax = theBAS2.LastVParameter();
 
     if( Precision::IsNegativeInfinite(VMin) ||
         Precision::IsPositiveInfinite(VMax))
@@ -473,7 +337,7 @@ static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
   }
 
   IntAna_QuadQuadGeo inter;
-  inter.Perform(aS1,aS2,Tolang,aTolSum, aHigh);
+  inter.Perform(aS1,aS2,Tolang,theTol, aHigh);
   if(inter.TypeInter() == IntAna_Ellipse)
   {
     const gp_Elips anEl = inter.Ellipse(1);
@@ -494,23 +358,19 @@ static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
 void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
                                 const TopoDS_Face& aF2)
 {
-  Standard_Boolean RestrictLine = Standard_False, hasCone = Standard_False;
-  
   if (myContext.IsNull()) {
     myContext=new IntTools_Context;
   }
 
   mySeqOfCurve.Clear();
-  myTolReached2d=0.;
-  myTolReached3d=0.;
   myIsDone = Standard_False;
   myNbrestr=0;//?
 
   myFace1=aF1;
   myFace2=aF2;
 
-  const BRepAdaptor_Surface aBAS1(myFace1, Standard_False);
-  const BRepAdaptor_Surface aBAS2(myFace2, Standard_False);
+  const BRepAdaptor_Surface& aBAS1 = myContext->SurfaceAdaptor(myFace1);
+  const BRepAdaptor_Surface& aBAS2 = myContext->SurfaceAdaptor(myFace2);
   GeomAbs_SurfaceType aType1=aBAS1.GetType();
   GeomAbs_SurfaceType aType2=aBAS2.GetType();
 
@@ -545,10 +405,12 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
   const Handle(Geom_Surface) S1=BRep_Tool::Surface(myFace1);
   const Handle(Geom_Surface) S2=BRep_Tool::Surface(myFace2);
 
-  const Standard_Real aTolF1=BRep_Tool::Tolerance(myFace1);
-  const Standard_Real aTolF2=BRep_Tool::Tolerance(myFace2);
+  Standard_Real aFuzz = myFuzzyValue / 2.;
+  myTolF1 = BRep_Tool::Tolerance(myFace1) + aFuzz;
+  myTolF2 = BRep_Tool::Tolerance(myFace2) + aFuzz;
+  myTol = myTolF1 + myTolF2;
 
-  Standard_Real TolArc = aTolF1 + aTolF2;
+  Standard_Real TolArc = myTol;
   Standard_Real TolTang = TolArc;
 
   const Standard_Boolean isFace1Quad = (aType1 == GeomAbs_Cylinder ||
@@ -562,43 +424,32 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
   if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane)  {
     Standard_Real umin, umax, vmin, vmax;
     //
-    BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
-    CorrectPlaneBoundaries(umin, umax, vmin, vmax);
+    myContext->UVBounds(myFace1, umin, umax, vmin, vmax);
     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
     //
-    BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
-    CorrectPlaneBoundaries(umin, umax, vmin, vmax);
+    myContext->UVBounds(myFace2, umin, umax, vmin, vmax);
     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
     //
     Standard_Real TolAng = 1.e-8;
     //
-    PerformPlanes(myHS1, myHS2, TolAng, TolTang, myApprox1, myApprox2, 
+    PerformPlanes(myHS1, myHS2,
+                  myTolF1, myTolF2, TolAng, TolTang,
+                  myApprox1, myApprox2,
                   mySeqOfCurve, myTangentFaces);
     //
     myIsDone = Standard_True;
-    
-    if(!myTangentFaces) {
+    //
+    if (!myTangentFaces) {
       const Standard_Integer NbLinPP = mySeqOfCurve.Length();
-      if(NbLinPP) {
-        Standard_Real aTolFMax;
-        myTolReached3d = 1.e-7;
-        aTolFMax=Max(aTolF1, aTolF2);
-        if (aTolFMax>myTolReached3d) {
-          myTolReached3d=aTolFMax;
-        }
-        //
-        myTolReached2d = myTolReached3d;
-
-        if (bReverse) {
-          Handle(Geom2d_Curve) aC2D1, aC2D2;
-          const Standard_Integer aNbLin = mySeqOfCurve.Length();
-          for (Standard_Integer i = 1; i <= aNbLin; ++i) {
-            IntTools_Curve& aIC=mySeqOfCurve(i);
-            aC2D1=aIC.FirstCurve2d();
-            aC2D2=aIC.SecondCurve2d();
-            aIC.SetFirstCurve2d(aC2D2);
-            aIC.SetSecondCurve2d(aC2D1);
-          }
+      if (NbLinPP && bReverse) {
+        Handle(Geom2d_Curve) aC2D1, aC2D2;
+        const Standard_Integer aNbLin = mySeqOfCurve.Length();
+        for (Standard_Integer i = 1; i <= aNbLin; ++i) {
+          IntTools_Curve& aIC = mySeqOfCurve(i);
+          aC2D1 = aIC.FirstCurve2d();
+          aC2D2 = aIC.SecondCurve2d();
+          aIC.SetFirstCurve2d(aC2D2);
+          aIC.SetSecondCurve2d(aC2D1);
         }
       }
     }
@@ -609,44 +460,34 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
   {
     Standard_Real umin, umax, vmin, vmax;
     // F1
-    BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax); 
+    myContext->UVBounds(myFace1, umin, umax, vmin, vmax); 
     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
     // F2
-    BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
-    CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
+    myContext->UVBounds(myFace2, umin, umax, vmin, vmax);
+    CorrectSurfaceBoundaries(myFace2, myTol * 2., umin, umax, vmin, vmax);
     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
-    //
-    if( aType2==GeomAbs_Cone ) { 
-      TolArc = 0.0001; 
-      hasCone = Standard_True; 
-    }
   }
   else if ((aType2==GeomAbs_Plane) && isFace1Quad)
   {
     Standard_Real umin, umax, vmin, vmax;
     //F1
-    BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
-    CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
+    myContext->UVBounds(myFace1, umin, umax, vmin, vmax);
+    CorrectSurfaceBoundaries(myFace1, myTol * 2., umin, umax, vmin, vmax);
     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
     // F2
-    BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
+    myContext->UVBounds(myFace2, umin, umax, vmin, vmax);
     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
-    //
-    if( aType1==GeomAbs_Cone ) {
-      TolArc = 0.0001; 
-      hasCone = Standard_True; 
-    }
   }
   else
   {
     Standard_Real umin, umax, vmin, vmax;
-    BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
-    CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
+    myContext->UVBounds(myFace1, umin, umax, vmin, vmax);
+    CorrectSurfaceBoundaries(myFace1, myTol * 2., umin, umax, vmin, vmax);
     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
-    BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
-    CorrectSurfaceBoundaries(myFace2, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
+    myContext->UVBounds(myFace2, umin, umax, vmin, vmax);
+    CorrectSurfaceBoundaries(myFace2, myTol * 2., umin, umax, vmin, vmax);
     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
   }
 
@@ -660,18 +501,10 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
 
   {
     const Standard_Real UVMaxStep = 0.001;
-    const Standard_Real Deflection = (hasCone) ? 0.085 : 0.1;
-  myIntersector.SetTolerances(TolArc, TolTang, UVMaxStep, Deflection); 
+    const Standard_Real Deflection = 0.1;
+    myIntersector.SetTolerances(TolArc, TolTang, UVMaxStep, Deflection); 
   }
   
-  if((myHS1->IsUClosed() && !myHS1->IsUPeriodic()) || 
-     (myHS1->IsVClosed() && !myHS1->IsVPeriodic()) ||
-     (myHS2->IsUClosed() && !myHS2->IsUPeriodic()) || 
-     (myHS2->IsVClosed() && !myHS2->IsVPeriodic()))
-  {
-    RestrictLine = Standard_True;
-  }
-  //
   if((aType1 != GeomAbs_BSplineSurface) &&
       (aType1 != GeomAbs_BezierSurface)  &&
      (aType1 != GeomAbs_OtherSurface)  &&
@@ -679,8 +512,6 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
       (aType2 != GeomAbs_BezierSurface)  &&
      (aType2 != GeomAbs_OtherSurface))
   {
-    RestrictLine = Standard_True;
-
     if ((aType1 == GeomAbs_Torus) ||
         (aType2 == GeomAbs_Torus))
     {
@@ -688,30 +519,31 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
     }
   }
 
-  //
-  if(!RestrictLine)
-  {
-    TopExp_Explorer aExp;
-    for(Standard_Integer i = 0; (!RestrictLine) && (i < 2); i++)
-    {
-      const TopoDS_Face& aF=(!i) ? myFace1 : myFace2;
-      aExp.Init(aF, TopAbs_EDGE);
-      for(; aExp.More(); aExp.Next())
+#ifdef INTTOOLS_FACEFACE_DEBUG
+    if(!myListOfPnts.IsEmpty()) {
+      char aBuff[10000];
+
+      Sprintf(aBuff,"bopcurves <face1 face2> -2d");
+      IntSurf_ListIteratorOfListOfPntOn2S IterLOP1(myListOfPnts);
+      for(;IterLOP1.More(); IterLOP1.Next())
       {
-        const TopoDS_Edge& aE=TopoDS::Edge(aExp.Current());
+        const IntSurf_PntOn2S& aPt = IterLOP1.Value();
+        Standard_Real u1, v1, u2, v2;
+        aPt.Parameters(u1, v1, u2, v2);
 
-        if(BRep_Tool::Degenerated(aE))
-        {
-          RestrictLine = Standard_True;
-          break;
-        }
+        Sprintf(aBuff, "%s -p %+10.20f %+10.20f %+10.20f %+10.20f", aBuff, u1, v1, u2, v2);
       }
+
+      cout << aBuff << endl;
     }
-  }
+#endif
 
-  const Standard_Boolean isGeomInt = isTreatAnalityc(aF1, aF2);
-  myIntersector.Perform(myHS1, dom1, myHS2, dom2, TolArc, TolTang, 
-                                  myListOfPnts, RestrictLine, isGeomInt);
+  const Standard_Boolean isGeomInt = isTreatAnalityc(aBAS1, aBAS2, myTol);
+  if (aF1.IsSame(aF2))
+    myIntersector.Perform(myHS1, dom1, TolArc, TolTang);
+  else
+    myIntersector.Perform(myHS1, dom1, myHS2, dom2, TolArc, TolTang, 
+                          myListOfPnts, isGeomInt);
 
   myIsDone = myIntersector.IsDone();
 
@@ -722,13 +554,9 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
       return;
     }
     //
-    if(RestrictLine) {
-      myListOfPnts.Clear(); // to use LineConstructor
-    }
-    //
-    const Standard_Integer aNbLin = myIntersector.NbLines();
-    for (Standard_Integer i=1; i <= aNbLin; ++i) {
-      MakeCurve(i, dom1, dom2);
+    const Standard_Integer aNbLinIntersector = myIntersector.NbLines();
+    for (Standard_Integer i=1; i <= aNbLinIntersector; ++i) {
+      MakeCurve(i, dom1, dom2, TolArc);
     }
     //
     ComputeTolReached3d();
@@ -736,8 +564,8 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
     if (bReverse) {
       Handle(Geom2d_Curve) aC2D1, aC2D2;
       //
-      const Standard_Integer aNbLin=mySeqOfCurve.Length();
-      for (Standard_Integer i=1; i<=aNbLin; ++i)
+      const Standard_Integer aNbLinSeqOfCurve =mySeqOfCurve.Length();
+      for (Standard_Integer i=1; i<=aNbLinSeqOfCurve; ++i)
       {
         IntTools_Curve& aIC=mySeqOfCurve(i);
         aC2D1=aIC.FirstCurve2d();
@@ -792,236 +620,76 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
 }
 
 //=======================================================================
-//function : ComputeTolerance
+//function :ComputeTolReached3d 
 //purpose  : 
 //=======================================================================
-Standard_Real IntTools_FaceFace::ComputeTolerance()
+void IntTools_FaceFace::ComputeTolReached3d()
 {
-  Standard_Integer i, j, aNbLin;
-  Standard_Real aFirst, aLast, aD, aDMax, aT;
-  Handle(Geom_Surface) aS1, aS2;
+  Standard_Integer i, j, aNbLin = mySeqOfCurve.Length();
+  if (!aNbLin) {
+    return;
+  }
   //
-  aDMax = 0;
-  aNbLin = mySeqOfCurve.Length();
+  // Minimal tangential tolerance for the curve
+  Standard_Real aTolFMax = Max(myTolF1, myTolF2);
   //
-  aS1 = myHS1->ChangeSurface().Surface();
-  aS2 = myHS2->ChangeSurface().Surface();
+  const Handle(Geom_Surface)& aS1 = myHS1->ChangeSurface().Surface();
+  const Handle(Geom_Surface)& aS2 = myHS2->ChangeSurface().Surface();
   //
   for (i = 1; i <= aNbLin; ++i)
   {
-    const IntTools_Curve& aIC = mySeqOfCurve(i);
+    IntTools_Curve& aIC = mySeqOfCurve(i);
     const Handle(Geom_Curve)& aC3D = aIC.Curve();
     if (aC3D.IsNull())
     {
       continue;
     }
     //
-    aFirst = aC3D->FirstParameter();
-    aLast  = aC3D->LastParameter();
+    Standard_Real aTolC = aIC.Tolerance();
+    Standard_Real aFirst = aC3D->FirstParameter();
+    Standard_Real aLast  = aC3D->LastParameter();
     //
+    // Compute the tolerance for the curve
     const Handle(Geom2d_Curve)& aC2D1 = aIC.FirstCurve2d();
     const Handle(Geom2d_Curve)& aC2D2 = aIC.SecondCurve2d();
     //
     for (j = 0; j < 2; ++j)
     {
       const Handle(Geom2d_Curve)& aC2D = !j ? aC2D1 : aC2D2;
-      const Handle(Geom_Surface)& aS = !j ? aS1 : aS2;
-      //
       if (!aC2D.IsNull())
       {
+        // Look for the maximal deviation between 3D and 2D curves
+        Standard_Real aD, aT;
+        const Handle(Geom_Surface)& aS = !j ? aS1 : aS2;
         if (IntTools_Tools::ComputeTolerance
             (aC3D, aC2D, aS, aFirst, aLast, aD, aT))
         {
-          if (aD > aDMax)
+          if (aD > aTolC)
           {
-            aDMax = aD;
+            aTolC = aD;
           }
         }
       }
-    }
-  }
-  //
-  return aDMax;
-}
-
-//=======================================================================
-//function :ComputeTolReached3d 
-//purpose  : 
-//=======================================================================
-  void IntTools_FaceFace::ComputeTolReached3d()
-{
-  Standard_Integer aNbLin, i;
-  GeomAbs_SurfaceType aType1, aType2;
-  //
-  aNbLin=myIntersector.NbLines();
-  if (!aNbLin) {
-    return;
-  }
-  //
-  aType1=myHS1->Surface().GetType();
-  aType2=myHS2->Surface().GetType();
-  //
-  if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
-    if (aNbLin==2){ 
-      Handle(IntPatch_Line) aIL1, aIL2;
-      IntPatch_IType aTL1, aTL2;
-      //
-      aIL1=myIntersector.Line(1);
-      aIL2=myIntersector.Line(2);
-      aTL1=aIL1->ArcType();
-      aTL2=aIL2->ArcType();
-      if (aTL1==IntPatch_Lin && aTL2==IntPatch_Lin) {
-        Standard_Real aD, aDTresh, dTol;
-        gp_Lin aL1, aL2;
-        //
-        dTol=1.e-8;
-        aDTresh=1.5e-6;
-        //
-        aL1=Handle(IntPatch_GLine)::DownCast(aIL1)->Line();
-        aL2=Handle(IntPatch_GLine)::DownCast(aIL2)->Line();
-        aD=aL1.Distance(aL2);
-        aD=0.5*aD;
-        if (aD<aDTresh) {
-          myTolReached3d=aD+dTol;
-        }
-        return;
-      }
-    }
-    //ZZ
-    if (aNbLin) {// Check the distances
-      Standard_Real aDMax;
-      //
-      aDMax = ComputeTolerance();
-      if (aDMax > 0.) {
-        myTolReached3d = aDMax;
-      }
-    }// if (aNbLin) 
-  }// if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
-  //
-  //904/G3 f
-  else if (aType1==GeomAbs_Plane && aType2==GeomAbs_Plane) {
-    Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
-    //
-    aTolTresh=1.e-7;
-    //
-    aTolF1 = BRep_Tool::Tolerance(myFace1);
-    aTolF2 = BRep_Tool::Tolerance(myFace2);
-    aTolFMax=Max(aTolF1, aTolF2);
-    //
-    if (aTolFMax>aTolTresh) {
-      myTolReached3d=aTolFMax;
-    }
-  }//if (aType1==GeomAbs_Plane && aType2==GeomAbs_Plane) {
-  //t
-  //IFV Bug OCC20297 
-  else if((aType1 == GeomAbs_Cylinder && aType2 == GeomAbs_Plane) ||
-          (aType2 == GeomAbs_Cylinder && aType1 == GeomAbs_Plane)) {
-    if(aNbLin == 1) {
-      const Handle(IntPatch_Line)& aIL1 = myIntersector.Line(1);
-      if(aIL1->ArcType() == IntPatch_Circle) {
-        gp_Circ aCir = Handle(IntPatch_GLine)::DownCast(aIL1)->Circle();
-        gp_XYZ aCirDir = aCir.Axis().Direction().XYZ();
-        gp_XYZ aPlDir;
-        gp_Pln aPln;
-        if(aType1 == GeomAbs_Plane) {
-          aPln = myHS1->Surface().Plane();
-        }
-        else {
-          aPln = myHS2->Surface().Plane();
-        }
-        aPlDir = aPln.Axis().Direction().XYZ();
-        Standard_Real cs = aCirDir*aPlDir;
-        if(cs < 0.) aPlDir.Reverse();
-        Standard_Real eps = 1.e-14;
-        if(!aPlDir.IsEqual(aCirDir, eps)) {
-          Standard_Integer aNbP = 11;
-          Standard_Real dt = 2.*M_PI / (aNbP - 1), t;
-          for(t = 0.; t < 2.*M_PI; t += dt) {
-            Standard_Real d = aPln.Distance(ElCLib::Value(t, aCir)); 
-            if(myTolReached3d < d) myTolReached3d = d;
-          }
-          myTolReached3d *= 1.1;
+      else
+      {
+        // Look for the maximal deviation between 3D curve and surface
+        const TopoDS_Face& aF = !j ? myFace1 : myFace2;
+        Standard_Real aD = FindMaxDistance(aC3D, aFirst, aLast, aF, myContext);
+        if (aD > aTolC)
+        {
+          aTolC = aD;
         }
-      } //aIL1->ArcType() == IntPatch_Circle
-    } //aNbLin == 1
-  } // aType1 == GeomAbs_Cylinder && aType2 == GeomAbs_Plane) 
-  //End IFV Bug OCC20297
-  //
-  else if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Torus) ||
-           (aType2==GeomAbs_Plane && aType1==GeomAbs_Torus)) {
-    aNbLin=mySeqOfCurve.Length();
-    if (aNbLin!=1) {
-      return;
-    }
-    //
-    Standard_Integer aNbP;
-    Standard_Real aT, aT1, aT2, dT, aUT, aVT, aUP, aVP;
-    Standard_Real aDP, aDT, aDmax;
-    gp_Pln aPln;
-    gp_Torus aTorus;
-    gp_Pnt aP, aPP, aPT;
-    //
-    const IntTools_Curve& aIC=mySeqOfCurve(1);
-    const Handle(Geom_Curve)& aC3D=aIC.Curve();
-    const Handle(Geom_BSplineCurve)& aBS=
-      Handle(Geom_BSplineCurve)::DownCast(aC3D);
-    if (aBS.IsNull()) {
-      return;
-    }
-    //
-    aT1=aBS->FirstParameter();
-    aT2=aBS->LastParameter();
-    //
-    aPln  =(aType1==GeomAbs_Plane) ? myHS1->Plane() : myHS2->Plane();
-    aTorus=(aType1==GeomAbs_Plane) ? myHS2->Torus() : myHS1->Torus();
-    //
-    aDmax=-1.;
-    aNbP=11;
-    dT=(aT2-aT1)/(aNbP-1);
-    for (i=0; i<aNbP; ++i) {
-      aT=aT1+i*dT;
-      if (i==aNbP-1) {
-        aT=aT2;
       }
-      //
-      aC3D->D0(aT, aP);
-      //
-      ElSLib::Parameters(aPln, aP, aUP, aVP);
-      aPP=ElSLib::Value(aUP, aVP, aPln);
-      aDP=aP.SquareDistance(aPP);
-      if (aDP>aDmax) {
-        aDmax=aDP;
-      }
-      //
-      ElSLib::Parameters(aTorus, aP, aUT, aVT);
-      aPT=ElSLib::Value(aUT, aVT, aTorus);
-      aDT=aP.SquareDistance(aPT);
-      if (aDT>aDmax) {
-        aDmax=aDT;
-      }
-    }
-    //
-    if (aDmax > myTolReached3d*myTolReached3d) {
-      myTolReached3d=sqrt(aDmax);
-      myTolReached3d=1.1*myTolReached3d;
     }
-  }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Torus) ||
-  //
-  else if ((aType1==GeomAbs_SurfaceOfRevolution && aType2==GeomAbs_Cylinder) ||
-           (aType2==GeomAbs_SurfaceOfRevolution && aType1==GeomAbs_Cylinder) ||
-           (aType1==GeomAbs_Plane && aType2==GeomAbs_Sphere) ||
-           (aType2==GeomAbs_Plane && aType1==GeomAbs_Sphere) ||
-           (aType1==GeomAbs_Plane && aType2==GeomAbs_SurfaceOfExtrusion) ||
-           (aType2==GeomAbs_Plane && aType1==GeomAbs_SurfaceOfExtrusion) ||
-           (aType1==GeomAbs_Plane && aType2==GeomAbs_BSplineSurface) ||
-           (aType2==GeomAbs_Plane && aType1==GeomAbs_BSplineSurface) ||
-           !myApprox) {
+    // Set the valid tolerance for the curve
+    aIC.SetTolerance(aTolC);
     //
-    Standard_Real aDMax;
-    //
-    aDMax = ComputeTolerance();
-    if (aDMax > myTolReached3d) {
-      myTolReached3d = aDMax;
+    // Set the tangential tolerance for the curve.
+    // Note, that, currently, computation of the tangential tolerance is
+    // implemented for the Plane/Plane case only.
+    // Thus, set the tangential tolerance equal to maximal tolerance of faces.
+    if (aIC.TangentialTolerance() < aTolFMax) {
+      aIC.SetTangentialTolerance(aTolFMax);
     }
   }
 }
@@ -1030,9 +698,10 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
 //function : MakeCurve
 //purpose  : 
 //=======================================================================
-  void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
-                                    const Handle(Adaptor3d_TopolTool)& dom1,
-                                    const Handle(Adaptor3d_TopolTool)& dom2) 
+void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
+                                  const Handle(Adaptor3d_TopolTool)& dom1,
+                                  const Handle(Adaptor3d_TopolTool)& dom2,
+                                  const Standard_Real theToler) 
 {
   Standard_Boolean bDone, rejectSurface, reApprox, bAvoidLineConstructor;
   Standard_Boolean ok, bPCurvesOk;
@@ -1043,8 +712,8 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
   IntPatch_IType typl;
   Handle(Geom_Curve) newc;
   //
-  const Standard_Real TOLCHECK   =0.0000001;
-  const Standard_Real TOLANGCHECK=0.1;
+  const Standard_Real TOLCHECK    = 1.e-7;
+  const Standard_Real TOLANGCHECK = 1.e-6;
   //
   rejectSurface = Standard_False;
   reApprox = Standard_False;
@@ -1059,25 +728,11 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
   typl = L->ArcType();
   //
   if(typl==IntPatch_Walking) {
-    Handle(IntPatch_Line) anewL;
-    //
-    const Handle(IntPatch_WLine)& aWLine=
-      Handle(IntPatch_WLine)::DownCast(L);
-    //DumpWLine(aWLine);
-
-    anewL = ComputePurgedWLine(aWLine);
-    if(anewL.IsNull()) {
+    Handle(IntPatch_WLine) aWLine (Handle(IntPatch_WLine)::DownCast(L));
+    if(aWLine.IsNull()) {
       return;
     }
-    L = anewL;
-    
-    //const Handle(IntPatch_WLine)& aWLineX = Handle(IntPatch_WLine)::DownCast(L);
-    //DumpWLine(aWLineX);
-
-    //
-    if(!myListOfPnts.IsEmpty()) {
-      bAvoidLineConstructor = Standard_True;
-    }
+    L = aWLine;
 
     Standard_Integer nbp = aWLine->NbPnts();
     const IntSurf_PntOn2S& p1 = aWLine->Point(1);
@@ -1093,6 +748,9 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
 
   typl=L->ArcType();
 
+  if(typl == IntPatch_Restriction)
+    bAvoidLineConstructor = Standard_True;
+
   //
   // Line Constructor
   if(!bAvoidLineConstructor) {
@@ -1138,11 +796,6 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
         new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
     }
     //
-    // myTolReached3d
-    if (typl == IntPatch_Lin) {
-      TolR3d (myFace1, myFace2, myTolReached3d);
-    }
-    //
     aNbParts=myLConstruct.NbParts();
     for (i=1; i<=aNbParts; i++) {
       Standard_Boolean bFNIt, bLPIt;
@@ -1159,46 +812,32 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
         Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
         aCurve.SetCurve(aCT3D);
         if (typl == IntPatch_Parabola) {
-          Standard_Real aTolF1, aTolF2, aTolBase;
-          
-          aTolF1 = BRep_Tool::Tolerance(myFace1);
-          aTolF2 = BRep_Tool::Tolerance(myFace2);
-          aTolBase=aTolF1+aTolF2;
-          myTolReached3d=IntTools_Tools::CurveTolerance(aCT3D, aTolBase);
+          Standard_Real aTolC = IntTools_Tools::CurveTolerance(aCT3D, myTol);
+          aCurve.SetTolerance(aTolC);
         }
         //
-        aCurve.SetCurve(new Geom_TrimmedCurve(newc, fprm, lprm));
         if(myApprox1) { 
           Handle (Geom2d_Curve) C2d;
           GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
                 myHS1->ChangeSurface().Surface(), newc, C2d);
-          if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
-            myTolReached2d=Tolpc;
-          }
-            //            
-            aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
-          }
-          else { 
-            Handle(Geom2d_BSplineCurve) H1;
-            //
-            aCurve.SetFirstCurve2d(H1);
-          }
+
+          if (C2d.IsNull())
+            continue;
+
+          aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d, fprm, lprm));
+        }
         //
         if(myApprox2) { 
           Handle (Geom2d_Curve) C2d;
           GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
                     myHS2->ChangeSurface().Surface(), newc, C2d);
-          if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
-            myTolReached2d=Tolpc;
-          }
-          //
-          aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
-          }
-        else { 
-          Handle(Geom2d_BSplineCurve) H1;
-          //
-          aCurve.SetSecondCurve2d(H1);
+
+          if (C2d.IsNull())
+            continue;
+
+          aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d, fprm, lprm));
         }
+        //
         mySeqOfCurve.Append(aCurve);
       } //if (!bFNIt && !bLPIt) {
       else {
@@ -1265,9 +904,6 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
         (Handle(IntPatch_GLine)::DownCast(L)->Ellipse());
     }
     //
-    // myTolReached3d
-    TolR3d (myFace1, myFace2, myTolReached3d);
-    //
     aNbParts=myLConstruct.NbParts();
     //
     Standard_Real aPeriod, aNul;
@@ -1291,14 +927,12 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
         else {
           gp_Pnt P1 = newc->Value(fprm);
           gp_Pnt P2 = newc->Value(aPeriod);
-          Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
-          aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
 
-          if(P1.Distance(P2) > aTolDist) {
+          if(P1.Distance(P2) > myTol) {
             Standard_Real anewpar = fprm;
 
             if(ParameterOutOfBoundary(fprm, newc, myFace1, myFace2, 
-                                      lprm, Standard_False, anewpar, myContext)) {
+                                      lprm, Standard_False, myTol, anewpar, myContext)) {
               fprm = anewpar;
             }
             aSeqFprm.Append(fprm);
@@ -1314,14 +948,12 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
         else {
           gp_Pnt P1 = newc->Value(aNul);
           gp_Pnt P2 = newc->Value(lprm);
-          Standard_Real aTolDist = BRep_Tool::Tolerance(myFace1) + BRep_Tool::Tolerance(myFace2);
-          aTolDist = (myTolReached3d > aTolDist) ? myTolReached3d : aTolDist;
 
-          if(P1.Distance(P2) > aTolDist) {
+          if(P1.Distance(P2) > myTol) {
             Standard_Real anewpar = lprm;
 
             if(ParameterOutOfBoundary(lprm, newc, myFace1, myFace2, 
-                                      fprm, Standard_True, anewpar, myContext)) {
+                                      fprm, Standard_True, myTol, anewpar, myContext)) {
               lprm = anewpar;
             }
             aSeqFprm.Append(aNul);
@@ -1335,7 +967,6 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
         aSeqLprm.Append(lprm);
       }
     }
-    
     //
     aNbParts=aSeqFprm.Length();
     for (i=1; i<=aNbParts; i++) {
@@ -1357,39 +988,17 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
             Handle (Geom2d_Curve) C2d;
             GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc, 
                         myHS1->ChangeSurface().Surface(), newc, C2d);
-            if(Tolpc>myTolReached2d || myTolReached2d==0) { 
-              myTolReached2d=Tolpc;
-            }
-            //
             aCurve.SetFirstCurve2d(C2d);
           }
-          else { //// 
-            Handle(Geom2d_BSplineCurve) H1;
-            aCurve.SetFirstCurve2d(H1);
-          }
-
 
           if(myApprox2) { 
             Handle (Geom2d_Curve) C2d;
             GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
                     myHS2->ChangeSurface().Surface(),newc,C2d);
-            if(Tolpc>myTolReached2d || myTolReached2d==0) { 
-              myTolReached2d=Tolpc;
-            }
-            //
             aCurve.SetSecondCurve2d(C2d);
           }
-          else { 
-            Handle(Geom2d_BSplineCurve) H1;
-            aCurve.SetSecondCurve2d(H1);
-          }
-        }
-        
-        else { 
-          Handle(Geom2d_BSplineCurve) H1;
-          aCurve.SetFirstCurve2d(H1);
-          aCurve.SetSecondCurve2d(H1);
         }
+        //
         mySeqOfCurve.Append(aCurve);
           //==============================================        
       } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
@@ -1410,31 +1019,16 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
               Handle (Geom2d_Curve) C2d;
               GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
                     myHS1->ChangeSurface().Surface(),newc,C2d);
-              if(Tolpc>myTolReached2d || myTolReached2d==0) { 
-                myTolReached2d=Tolpc;
-              }
-              //
               aCurve.SetFirstCurve2d(C2d);
             }
-            else { //// 
-              Handle(Geom2d_BSplineCurve) H1;
-              aCurve.SetFirstCurve2d(H1);
-            }
 
             if(myApprox2) { 
               Handle (Geom2d_Curve) C2d;
               GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
                     myHS2->ChangeSurface().Surface(),newc,C2d);
-              if(Tolpc>myTolReached2d || myTolReached2d==0) { 
-                myTolReached2d=Tolpc;
-              }
-              //
               aCurve.SetSecondCurve2d(C2d);
             }
-            else { 
-              Handle(Geom2d_BSplineCurve) H1;
-              aCurve.SetSecondCurve2d(H1);
-            }
+            //
             mySeqOfCurve.Append(aCurve);
             break;
           }
@@ -1463,40 +1057,16 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
                 Handle (Geom2d_Curve) C2d;
                 GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc, 
                         myHS1->ChangeSurface().Surface(), newc, C2d);
-                if(Tolpc>myTolReached2d || myTolReached2d==0) { 
-                  myTolReached2d=Tolpc;
-                }
-                // 
                 aCurve.SetFirstCurve2d(C2d);
               }
-              else { 
-                Handle(Geom2d_BSplineCurve) H1;
-                aCurve.SetFirstCurve2d(H1);
-              }
-                
+
               if(myApprox2) { 
                 Handle (Geom2d_Curve) C2d;
                 GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
                         myHS2->ChangeSurface().Surface(), newc, C2d);
-                if(Tolpc>myTolReached2d || myTolReached2d==0) { 
-                  myTolReached2d=Tolpc;
-                }
-                //                 
                 aCurve.SetSecondCurve2d(C2d);
               }
-                
-              else { 
-                Handle(Geom2d_BSplineCurve) H1;
-                aCurve.SetSecondCurve2d(H1);
-              }
             }//  end of if (typl == IntPatch_Circle || typl == IntPatch_Ellipse)
-             
-            else { 
-              Handle(Geom2d_BSplineCurve) H1;
-              //         
-              aCurve.SetFirstCurve2d(H1);
-              aCurve.SetSecondCurve2d(H1);
-            }
             //==============================================        
             //
             mySeqOfCurve.Append(aCurve);
@@ -1509,188 +1079,18 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
     }// IntPatch_Circle: IntPatch_Ellipse:
     break;
     
-  case IntPatch_Analytic: {
-    IntSurf_Quadric quad1,quad2;
-    GeomAbs_SurfaceType typs = myHS1->Surface().GetType();
-    
-    switch (typs) {
-      case GeomAbs_Plane:
-        quad1.SetValue(myHS1->Surface().Plane());
-        break;
-      case GeomAbs_Cylinder:
-        quad1.SetValue(myHS1->Surface().Cylinder());
-        break;
-      case GeomAbs_Cone:
-        quad1.SetValue(myHS1->Surface().Cone());
-        break;
-      case GeomAbs_Sphere:
-        quad1.SetValue(myHS1->Surface().Sphere());
-        break;
-    case GeomAbs_Torus:
-      quad1.SetValue(myHS1->Surface().Torus());
-      break;
-      default:
-        Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 1");
-      }
-      
-    typs = myHS2->Surface().GetType();
-    
-    switch (typs) {
-      case GeomAbs_Plane:
-        quad2.SetValue(myHS2->Surface().Plane());
-        break;
-      case GeomAbs_Cylinder:
-        quad2.SetValue(myHS2->Surface().Cylinder());
-        break;
-      case GeomAbs_Cone:
-        quad2.SetValue(myHS2->Surface().Cone());
-        break;
-      case GeomAbs_Sphere:
-        quad2.SetValue(myHS2->Surface().Sphere());
-        break;
-    case GeomAbs_Torus:
-      quad2.SetValue(myHS2->Surface().Torus());
-      break;
-      default:
-        Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 2");
-      }
-    //
-    //=========
-    IntPatch_ALineToWLine convert (quad1, quad2);
-      
-    if (!myApprox) {
-      aNbParts=myLConstruct.NbParts();
-      for (i=1; i<=aNbParts; i++) {
-        myLConstruct.Part(i, fprm, lprm);
-        Handle(IntPatch_WLine) WL = 
-          convert.MakeWLine(Handle(IntPatch_ALine)::DownCast(L), fprm, lprm);
-        //
-        Handle(Geom2d_BSplineCurve) H1;
-        Handle(Geom2d_BSplineCurve) H2;
-
-        if(myApprox1) {
-          H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
-        }
-        
-        if(myApprox2) {
-          H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
-        }
-        //          
-        mySeqOfCurve.Append(IntTools_Curve(MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
-      }
-    } // if (!myApprox)
-
-    else { // myApprox=TRUE
-      GeomInt_WLApprox theapp3d;
-      // 
-      Standard_Real tol2d = myTolApprox;
-      //         
-      theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
-      
-      aNbParts=myLConstruct.NbParts();
-      for (i=1; i<=aNbParts; i++) {
-        myLConstruct.Part(i, fprm, lprm);
-        Handle(IntPatch_WLine) WL = 
-          convert.MakeWLine(Handle(IntPatch_ALine):: DownCast(L),fprm,lprm);
-
-        theapp3d.Perform(myHS1,myHS2,WL,Standard_True,myApprox1,myApprox2, 1, WL->NbPnts());
-
-        if (!theapp3d.IsDone()) {
-          //
-          Handle(Geom2d_BSplineCurve) H1;
-          Handle(Geom2d_BSplineCurve) H2;
-
-          if(myApprox1) {
-            H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
-          }
-          
-          if(myApprox2) {
-            H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
-          }
-          //          
-          mySeqOfCurve.Append(IntTools_Curve(MakeBSpline(WL,1,WL->NbPnts()), H1, H2));
-        }
-
-        else {
-          if(myApprox1 || myApprox2) { 
-            if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0) { 
-              myTolReached2d = theapp3d.TolReached2d();
-            }
-          }
-          
-          if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0) { 
-            myTolReached3d = theapp3d.TolReached3d();
-          }
-
-          Standard_Integer aNbMultiCurves, nbpoles;
-          aNbMultiCurves=theapp3d.NbMultiCurves();
-          for (j=1; j<=aNbMultiCurves; j++) {
-            const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
-            nbpoles = mbspc.NbPoles();
-            
-            TColgp_Array1OfPnt tpoles(1, nbpoles);
-            mbspc.Curve(1, tpoles);
-            Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
-                                                               mbspc.Knots(),
-                                                               mbspc.Multiplicities(),
-                                                               mbspc.Degree());
-            
-            GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
-            Check.FixTangent(Standard_True,Standard_True);
-            // 
-            IntTools_Curve aCurve;
-            aCurve.SetCurve(BS);
-            
-            if(myApprox1) { 
-              TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
-              mbspc.Curve(2,tpoles2d);
-              Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
-                                                                      mbspc.Knots(),
-                                                                      mbspc.Multiplicities(),
-                                                                      mbspc.Degree());
-
-              GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
-              newCheck.FixTangent(Standard_True,Standard_True);
-              //                 
-              aCurve.SetFirstCurve2d(BS2);
-            }
-            else {
-              Handle(Geom2d_BSplineCurve) H1;
-              aCurve.SetFirstCurve2d(H1);
-            }
-            
-            if(myApprox2) { 
-              TColgp_Array1OfPnt2d tpoles2d(1, nbpoles);
-              Standard_Integer TwoOrThree;
-              TwoOrThree=myApprox1 ? 3 : 2;
-              mbspc.Curve(TwoOrThree, tpoles2d);
-              Handle(Geom2d_BSplineCurve) BS2 =new Geom2d_BSplineCurve(tpoles2d,
-                                                                       mbspc.Knots(),
-                                                                       mbspc.Multiplicities(),
-                                                                       mbspc.Degree());
-                
-              GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
-              newCheck.FixTangent(Standard_True,Standard_True);
-              //         
-              aCurve.SetSecondCurve2d(BS2);
-            }
-            else { 
-              Handle(Geom2d_BSplineCurve) H2;
-              aCurve.SetSecondCurve2d(H2);
-            }
-            // 
-            mySeqOfCurve.Append(aCurve);
-
-          }// for (j=1; j<=aNbMultiCurves; j++) {
-        }// else from if (!theapp3d.IsDone())
-      }// for (i=1; i<=aNbParts; i++) {
-    }// else { // myApprox=TRUE
-  }// case IntPatch_Analytic:
+  case IntPatch_Analytic:
+    //This case was processed earlier (in IntPatch_Intersection)
     break;
 
   case IntPatch_Walking:{
     Handle(IntPatch_WLine) WL = 
       Handle(IntPatch_WLine)::DownCast(L);
+
+#ifdef INTTOOLS_FACEFACE_DEBUG
+    WL->Dump(0);
+#endif
+
     //
     Standard_Integer ifprm, ilprm;
     //
@@ -1714,14 +1114,14 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
         }
         //
         if(myApprox1) {
-          H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
+          H1 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
         }
         //
         if(myApprox2) {
-          H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
+          H2 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
         }
         //           
-        aBSp=MakeBSpline(WL, ifprm, ilprm);
+        aBSp=GeomInt_IntSS::MakeBSpline(WL, ifprm, ilprm);
         IntTools_Curve aIC(aBSp, H1, H2);
         mySeqOfCurve.Append(aIC);
       }// for (i=1; i<=aNbParts; ++i) {
@@ -1764,39 +1164,41 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
       }
         
       if(myHS1 == myHS2) { 
-        theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
+        theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30, Standard_False, aParType);
         rejectSurface = Standard_True;
       }
       else { 
         if(reApprox && !rejectSurface)
-          theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
+          theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30, Standard_False, aParType);
         else {
           Standard_Integer iDegMax, iDegMin, iNbIter;
           //
           ApproxParameters(myHS1, myHS2, iDegMin, iDegMax, iNbIter);
-          theapp3d.SetParameters(myTolApprox, tol2d, iDegMin, iDegMax, iNbIter, Standard_True, aParType);
+          theapp3d.SetParameters(myTolApprox, tol2d, iDegMin, iDegMax,
+                                                  iNbIter, 30, Standard_True, aParType);
         }
       }
       //
       Standard_Real aReachedTol = Precision::Confusion();
-      bIsDecomposited=DecompositionOfWLine(WL,
-                                           myHS1, 
-                                           myHS2, 
-                                           myFace1, 
-                                           myFace2, 
-                                           myLConstruct, 
-                                           bAvoidLineConstructor, 
-                                           aSeqOfL, 
-                                           aReachedTol,
-                                           myContext);
-      if ( bIsDecomposited && ( myTolReached3d < aReachedTol ) ) {
-        myTolReached3d = aReachedTol;
-      }
+      bIsDecomposited = IntTools_WLineTool::
+        DecompositionOfWLine(WL,
+                             myHS1, 
+                             myHS2, 
+                             myFace1, 
+                             myFace2, 
+                             myLConstruct, 
+                             bAvoidLineConstructor, 
+                             myTol,
+                             aSeqOfL, 
+                             aReachedTol,
+                             myContext);
       //
       aNbSeqOfL=aSeqOfL.Length();
       //
+      Standard_Real aTolC = 0.;
       if (bIsDecomposited) {
         nbiter=aNbSeqOfL;
+        aTolC = aReachedTol;
       }
       else {
         nbiter=1;
@@ -1839,15 +1241,17 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
                 (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
              
-              theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True, aParType);
+              theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30,
+                                                                Standard_True, aParType);
               
               Standard_Boolean bUseSurfaces;
-              bUseSurfaces=NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
+              bUseSurfaces = IntTools_WLineTool::NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
               if (bUseSurfaces) {
                 // ######
                 rejectSurface = Standard_True;
                 // ######
-                theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
+                theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30,
+                                                                Standard_False, aParType);
               }
             }
           }
@@ -1859,40 +1263,29 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
           Handle(Geom2d_BSplineCurve) H1;
           Handle(Geom2d_BSplineCurve) H2;
           //           
-          Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
+          Handle(Geom_Curve) aBSp=GeomInt_IntSS::MakeBSpline(WL,ifprm, ilprm);
           //
           if(myApprox1) {
-            H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
+            H1 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
           }
           //
           if(myApprox2) {
-            H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
+            H2 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
           }
           //           
           IntTools_Curve aIC(aBSp, H1, H2);
           mySeqOfCurve.Append(aIC);
         }
-        
         else {
-          if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) { 
-            if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) { 
-              myTolReached2d = theapp3d.TolReached2d();
-            }
-          }
-          if(typs1==GeomAbs_Plane || typs2==GeomAbs_Plane) { 
-            myTolReached3d = myTolReached2d;
+          if (typs1 == GeomAbs_Plane || typs2 == GeomAbs_Plane) {
             //
-            if (typs1==GeomAbs_Torus || typs2==GeomAbs_Torus) {
-              if (myTolReached3d<1.e-6) {
-                myTolReached3d = theapp3d.TolReached3d();
-                myTolReached3d=1.e-6;
+            if (typs1 == GeomAbs_Torus || typs2 == GeomAbs_Torus) {
+              if (aTolC < 1.e-6) {
+                aTolC = 1.e-6;
               }
             }
           }
-          else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) { 
-            myTolReached3d = theapp3d.TolReached3d();
-          }
-          
+          //
           Standard_Integer aNbMultiCurves, nbpoles;
           aNbMultiCurves=theapp3d.NbMultiCurves(); 
           for (j=1; j<=aNbMultiCurves; j++) {
@@ -1945,10 +1338,6 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
                 // ############################################
                 aCurve.SetFirstCurve2d(BS1);
               }
-              else {
-                Handle(Geom2d_BSplineCurve) H1;
-                aCurve.SetFirstCurve2d(H1);
-              }
 
               if(myApprox2) { 
                 mbspc.Curve(2, tpoles2d);
@@ -1972,16 +1361,15 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
                 // 
                 aCurve.SetSecondCurve2d(BS2);
               }
-              else { 
-                Handle(Geom2d_BSplineCurve) H2;
-                //                 
-                  aCurve.SetSecondCurve2d(H2);
-              }
+              //
+              aCurve.SetTolerance(aTolC);
               //
               mySeqOfCurve.Append(aCurve);
+
             }//if(typs1 == GeomAbs_Plane) {
             
-            else if(typs2 == GeomAbs_Plane) { 
+            else if(typs2 == GeomAbs_Plane)
+            {
               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
               nbpoles = mbspc.NbPoles();
               
@@ -2008,6 +1396,7 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
               //         
               IntTools_Curve aCurve;
               aCurve.SetCurve(BS);
+              aCurve.SetTolerance(aTolC);
 
               if(myApprox2) {
                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
@@ -2026,14 +1415,10 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
                   }
                 }
                 // ###########################################
-                bPCurvesOk = CheckPCurve(BS1, myFace2);
+                bPCurvesOk = CheckPCurve(BS1, myFace2, myContext);
                 aCurve.SetSecondCurve2d(BS1);
               }
-              else {
-                Handle(Geom2d_BSplineCurve) H2;
-                aCurve.SetSecondCurve2d(H2);
-              }
-              
+
               if(myApprox1) { 
                 mbspc.Curve(1,tpoles2d);
                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
@@ -2052,14 +1437,9 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
                   }
                 }
                 // ###########################################
-                bPCurvesOk = bPCurvesOk && CheckPCurve(BS2, myFace1);
+                bPCurvesOk = bPCurvesOk && CheckPCurve(BS2, myFace1, myContext);
                 aCurve.SetFirstCurve2d(BS2);
               }
-              else { 
-                Handle(Geom2d_BSplineCurve) H1;
-                //                 
-                aCurve.SetFirstCurve2d(H1);
-              }
               //
               //if points of the pcurves are out of the faces bounds
               //create 3d and 2d curves without approximation
@@ -2067,22 +1447,22 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
                 Handle(Geom2d_BSplineCurve) H1, H2;
                 bPCurvesOk = Standard_True;
                 //           
-                Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
+                Handle(Geom_Curve) aBSp=GeomInt_IntSS::MakeBSpline(WL,ifprm, ilprm);
                 
                 if(myApprox1) {
-                  H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
-                  bPCurvesOk = CheckPCurve(H1, myFace1);
+                  H1 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
+                  bPCurvesOk = CheckPCurve(H1, myFace1, myContext);
                 }
                 
                 if(myApprox2) {
-                  H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
-                  bPCurvesOk = bPCurvesOk && CheckPCurve(H2, myFace2);
+                  H2 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
+                  bPCurvesOk = bPCurvesOk && CheckPCurve(H2, myFace2, myContext);
                 }
                 //
                 //if pcurves created without approximation are out of the 
                 //faces bounds, use approximated 3d and 2d curves
                 if (bPCurvesOk) {
-                  IntTools_Curve aIC(aBSp, H1, H2);
+                  IntTools_Curve aIC(aBSp, H1, H2, aTolC);
                   mySeqOfCurve.Append(aIC);
                 } else {
                   mySeqOfCurve.Append(aCurve);
@@ -2090,13 +1470,13 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
               } else {
                 mySeqOfCurve.Append(aCurve);
               }
+
             }// else if(typs2 == GeomAbs_Plane)
             //
             else { //typs2 != GeomAbs_Plane && typs1 != GeomAbs_Plane
               Standard_Boolean bIsValid1, bIsValid2;
               Handle(Geom_BSplineCurve) BS;
-              Handle(Geom2d_BSplineCurve) aH2D;        
-              IntTools_Curve aCurve; 
+              IntTools_Curve aCurve;
               //
               bIsValid1=Standard_True;
               bIsValid2=Standard_True;
@@ -2111,10 +1491,9 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
                                                                  mbspc.Degree());
               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
               Check.FixTangent(Standard_True,Standard_True);
-              //                 
+              //
               aCurve.SetCurve(BS);
-              aCurve.SetFirstCurve2d(aH2D);
-              aCurve.SetSecondCurve2d(aH2D);
+              aCurve.SetTolerance(aTolC);
               //
               if(myApprox1) { 
                 if(anApprox1) {
@@ -2123,14 +1502,14 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
                   mbspc.Curve(2,tpoles2d);
                   //
                   BS1=new Geom2d_BSplineCurve(tpoles2d,
-                                                                        mbspc.Knots(),
-                                                                        mbspc.Multiplicities(),
-                                                                        mbspc.Degree());
+                                              mbspc.Knots(),
+                                              mbspc.Multiplicities(),
+                                              mbspc.Degree());
                   GeomLib_Check2dBSplineCurve newCheck(BS1,TOLCHECK,TOLANGCHECK);
                   newCheck.FixTangent(Standard_True,Standard_True);
                   //         
                   if (!reApprox) {
-                    bIsValid1=CheckPCurve(BS1, myFace1);
+                    bIsValid1=CheckPCurve(BS1, myFace1, myContext);
                   }
                   //
                   aCurve.SetFirstCurve2d(BS1);
@@ -2162,7 +1541,7 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
                   newCheck.FixTangent(Standard_True,Standard_True);
                 //                 
                   if (!reApprox) {
-                    bIsValid2=CheckPCurve(BS2, myFace2);        
+                    bIsValid2=CheckPCurve(BS2, myFace2, myContext);
                   }
                   aCurve.SetSecondCurve2d(BS2);
                 }
@@ -2197,32 +1576,13 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
     
   case IntPatch_Restriction: 
     {
-      GeomAbs_SurfaceType typS1 = myHS1->Surface().GetType();
-      GeomAbs_SurfaceType typS2 = myHS2->Surface().GetType();
-      Standard_Boolean isAnalS1 = Standard_False;
-      switch (typS1)
-      {
-      case GeomAbs_Plane:
-      case GeomAbs_Cylinder:
-      case GeomAbs_Sphere:
-      case GeomAbs_Cone: 
-      case GeomAbs_Torus: isAnalS1 = Standard_True; break;
-      default: break;
-      }
-
-      Standard_Integer isAnalS2 = Standard_False;
-      switch (typS2)
-      {
-      case GeomAbs_Plane:
-      case GeomAbs_Cylinder:
-      case GeomAbs_Sphere:
-      case GeomAbs_Cone: 
-      case GeomAbs_Torus: isAnalS2 = Standard_True; break;
-      default: break;
-      }
-
       Handle(IntPatch_RLine) RL = 
         Handle(IntPatch_RLine)::DownCast(L);
+
+#ifdef INTTOOLS_FACEFACE_DEBUG
+    RL->Dump(0);
+#endif
+
       Handle(Geom_Curve) aC3d;
       Handle(Geom2d_Curve) aC2d1, aC2d2;
       Standard_Real aTolReached;
@@ -2257,16 +1617,27 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
       GeomInt_IntSS::
         TrimILineOnSurfBoundaries(aC2d1, aC2d2, aBox1, aBox2, anArrayOfParameters);
 
+      //Intersect with true boundaries. After that, enlarge bounding-boxes in order to 
+      //correct definition, if point on curve is inscribed in the box.
+      aBox1.Enlarge(theToler);
+      aBox2.Enlarge(theToler);
+
       const Standard_Integer aNbIntersSolutionsm1 = anArrayOfParameters.Length() - 1;
 
       //Trim RLine found.
       for(Standard_Integer anInd = 0; anInd < aNbIntersSolutionsm1; anInd++)
       {
-        const Standard_Real aParF = anArrayOfParameters(anInd),
-                            aParL = anArrayOfParameters(anInd+1);
+        Standard_Real &aParF = anArrayOfParameters(anInd),
+                      &aParL = anArrayOfParameters(anInd+1);
 
         if((aParL - aParF) <= Precision::PConfusion())
+        {
+          //In order to more precise extending to the boundaries of source curves.
+          if(anInd < aNbIntersSolutionsm1-1)
+            aParL = aParF;
+
           continue;
+        }
 
         const Standard_Real aPar = 0.5*(aParF + aParL);
         gp_Pnt2d aPt;
@@ -2340,7 +1711,7 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
     quad1.SetValue(HS1->Surface().Torus());
     break;
   default:
-    Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
+    throw Standard_ConstructionError("GeomInt_IntSS::MakeCurve");
   }
   
   typs = HS2->Surface().GetType();
@@ -2361,7 +1732,7 @@ Standard_Real IntTools_FaceFace::ComputeTolerance()
     quad2.SetValue(HS2->Surface().Torus());
     break;
   default:
-    Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
+    throw Standard_ConstructionError("GeomInt_IntSS::MakeCurve");
   }
 
   quad1.Parameters(Ptref,U1,V1);
@@ -2390,37 +1761,7 @@ Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
   return
     new Geom_BSplineCurve(poles,knots,mults,1);
 }
-//
-
-//=======================================================================
-//function : MakeBSpline2d
-//purpose  : 
-//=======================================================================
-Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
-                                          const Standard_Integer ideb,
-                                          const Standard_Integer ifin,
-                                          const Standard_Boolean onFirst)
-{
-  Standard_Integer i, nbpnt = ifin-ideb+1;
-  TColgp_Array1OfPnt2d poles(1,nbpnt);
-  TColStd_Array1OfReal knots(1,nbpnt);
-  TColStd_Array1OfInteger mults(1,nbpnt);
-  Standard_Integer ipidebm1;
-
-  for(i = 1, ipidebm1 = i+ideb-1; i <= nbpnt; ipidebm1++, i++) {
-      Standard_Real U, V;
-      if(onFirst)
-        theWLine->Point(ipidebm1).ParametersOnS1(U, V);
-      else
-        theWLine->Point(ipidebm1).ParametersOnS2(U, V);
-      poles(i).SetCoord(U, V);
-      mults(i) = 1;
-      knots(i) = i-1;
-    }
-    mults(1) = mults(nbpnt) = 2;
 
-  return new Geom2d_BSplineCurve(poles,knots,mults,1);
-}
 //=======================================================================
 //function : PrepareLines3D
 //purpose  : 
@@ -2583,59 +1924,6 @@ Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine
     }
   }
   //
-  {
-    Standard_Integer aNbP;
-    Standard_Real aXP, dXfact, aXmid, aX1, aX2, aTolPA;
-    //
-    aTolPA=Precision::Angular();
-    // U
-    if (isuperiodic) {
-      aXP=anAdaptorSurface.UPeriod();
-      dXfact=theumax-theumin;
-      if (dXfact-aTolPA>aXP) {
-        aXmid=0.5*(theumax+theumin);
-        aNbP=RealToInt(aXmid/aXP);
-        if (aXmid<0.) {
-          aNbP=aNbP-1;
-        }
-        aX1=aNbP*aXP;
-        if (theumin>aTolPA) {
-          aX1=theumin+aNbP*aXP;
-        }
-        aX2=aX1+aXP;
-        if (theumin<aX1) {
-          theumin=aX1;
-        }
-        if (theumax>aX2) {
-          theumax=aX2;
-        }
-      }
-    }
-    // V
-    if (isvperiodic) {
-      aXP=anAdaptorSurface.VPeriod();
-      dXfact=thevmax-thevmin;
-      if (dXfact-aTolPA>aXP) {
-        aXmid=0.5*(thevmax+thevmin);
-        aNbP=RealToInt(aXmid/aXP);
-        if (aXmid<0.) {
-          aNbP=aNbP-1;
-        }
-        aX1=aNbP*aXP;
-        if (thevmin>aTolPA) {
-          aX1=thevmin+aNbP*aXP;
-        }
-        aX2=aX1+aXP;
-        if (thevmin<aX1) {
-          thevmin=aX1;
-        }
-        if (thevmax>aX2) {
-          thevmax=aX2;
-        }
-      }
-    }
-  }
-  //
   if(isuperiodic || isvperiodic) {
     Standard_Boolean correct = Standard_False;
     Standard_Boolean correctU = Standard_False;
@@ -2690,16 +1978,13 @@ Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine
       aBox.Get(umin, vmin, umax, vmax);
 
       if(isuperiodic && correctU) {
-        
         if(theumin < umin)
           theumin = umin;
-        
         if(theumax > umax) {
           theumax = umax;
         }
       }
       if(isvperiodic && correctV) {
-        
         if(thevmin < vmin)
           thevmin = vmin;
         if(thevmax > vmax)
@@ -2708,1558 +1993,80 @@ Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine
     }
   }
 }
-//
-//
-// The block is dedicated to determine whether WLine [ifprm, ilprm]
-// crosses the degenerated zone on each given surface or not.
-// If Yes -> We will not use info about surfaces during approximation
-// because inside degenerated zone of the surface the approx. algo.
-// uses wrong values of normal, etc., and resulting curve will have
-// oscillations that we would not like to have. 
-
 
-static
-  Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
-                                     const Handle(Geom_Surface)& aS,
-                                     const Standard_Integer iDir);
-static
-  Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
-                                            const TopoDS_Face& aF1,
-                                            const TopoDS_Face& aF2);
-//=======================================================================
-//function :  NotUseSurfacesForApprox
-//purpose  : 
-//=======================================================================
-Standard_Boolean NotUseSurfacesForApprox(const TopoDS_Face& aF1,
-                                         const TopoDS_Face& aF2,
-                                         const Handle(IntPatch_WLine)& WL,
-                                         const Standard_Integer ifprm,
-                                         const Standard_Integer ilprm)
+// ------------------------------------------------------------------------------------------------
+// static function: ParameterOutOfBoundary
+// purpose:         Computes a new parameter for given curve. The corresponding 2d points 
+//                  does not lay on any boundary of given faces
+// ------------------------------------------------------------------------------------------------
+Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter, 
+                                        const Handle(Geom_Curve)& theCurve, 
+                                        const TopoDS_Face&        theFace1, 
+                                        const TopoDS_Face&        theFace2,
+                                        const Standard_Real       theOtherParameter,
+                                        const Standard_Boolean    bIncreasePar,
+                                        const Standard_Real       theTol,
+                                        Standard_Real&            theNewParameter,
+                                        const Handle(IntTools_Context)& aContext)
 {
-  Standard_Boolean bPInDZ;
+  Standard_Boolean bIsComputed = Standard_False;
+  theNewParameter = theParameter;
 
-  Handle(IntSurf_LineOn2S) aLineOn2S=WL->Curve();
-  
-  const IntSurf_PntOn2S& aP2Sfprm=aLineOn2S->Value(ifprm);
-  bPInDZ=IsPointInDegeneratedZone(aP2Sfprm, aF1, aF2);
-  if (bPInDZ) {
-    return bPInDZ;
-  }
+  Standard_Real acurpar = theParameter;
+  TopAbs_State aState = TopAbs_ON;
+  Standard_Integer iter = 0;
+  Standard_Real asumtol = theTol;
+  Standard_Real adelta = asumtol * 0.1;
+  adelta = (adelta < Precision::Confusion()) ? Precision::Confusion() : adelta;
+  Handle(Geom_Surface) aSurf1 = BRep_Tool::Surface(theFace1);
+  Handle(Geom_Surface) aSurf2 = BRep_Tool::Surface(theFace2);
 
-  const IntSurf_PntOn2S& aP2Slprm=aLineOn2S->Value(ilprm);
-  bPInDZ=IsPointInDegeneratedZone(aP2Slprm, aF1, aF2);
-  
-  return bPInDZ;
-}
-//=======================================================================
-//function : IsPointInDegeneratedZone
-//purpose  : 
-//=======================================================================
-Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
-                                          const TopoDS_Face& aF1,
-                                          const TopoDS_Face& aF2)
-                                          
-{
-  Standard_Boolean bFlag=Standard_True;
-  Standard_Real US11, US12, VS11, VS12, US21, US22, VS21, VS22;
-  Standard_Real U1, V1, U2, V2, aDelta, aD;
-  gp_Pnt2d aP2d;
-
-  Handle(Geom_Surface)aS1 = BRep_Tool::Surface(aF1);
-  aS1->Bounds(US11, US12, VS11, VS12);
-  GeomAdaptor_Surface aGAS1(aS1);
-
-  Handle(Geom_Surface)aS2 = BRep_Tool::Surface(aF2);
-  aS1->Bounds(US21, US22, VS21, VS22);
-  GeomAdaptor_Surface aGAS2(aS2);
-  //
-  //const gp_Pnt& aP=aP2S.Value();
-  aP2S.Parameters(U1, V1, U2, V2);
-  //
-  aDelta=1.e-7;
-  // Check on Surf 1
-  aD=aGAS1.UResolution(aDelta);
-  aP2d.SetCoord(U1, V1);
-  if (fabs(U1-US11) < aD) {
-    bFlag=IsDegeneratedZone(aP2d, aS1, 1);
-    if (bFlag) {
-      return bFlag;
+  Standard_Real u1, u2, v1, v2;
+
+  GeomAPI_ProjectPointOnSurf aPrj1;
+  aSurf1->Bounds(u1, u2, v1, v2);
+  aPrj1.Init(aSurf1, u1, u2, v1, v2);
+
+  GeomAPI_ProjectPointOnSurf aPrj2;
+  aSurf2->Bounds(u1, u2, v1, v2);
+  aPrj2.Init(aSurf2, u1, u2, v1, v2);
+
+  while(aState == TopAbs_ON) {
+    if(bIncreasePar)
+      acurpar += adelta;
+    else
+      acurpar -= adelta;
+    gp_Pnt aPCurrent = theCurve->Value(acurpar);
+    aPrj1.Perform(aPCurrent);
+    Standard_Real U=0., V=0.;
+
+    if(aPrj1.IsDone()) {
+      aPrj1.LowerDistanceParameters(U, V);
+      aState = aContext->StatePointFace(theFace1, gp_Pnt2d(U, V));
     }
-  }
-  if (fabs(U1-US12) < aD) {
-    bFlag=IsDegeneratedZone(aP2d, aS1, 1);
-    if (bFlag) {
-      return bFlag;
+
+    if(aState != TopAbs_ON) {
+      aPrj2.Perform(aPCurrent);
+                
+      if(aPrj2.IsDone()) {
+        aPrj2.LowerDistanceParameters(U, V);
+        aState = aContext->StatePointFace(theFace2, gp_Pnt2d(U, V));
+      }
     }
-  }
-  aD=aGAS1.VResolution(aDelta);
-  if (fabs(V1-VS11) < aDelta) {
-    bFlag=IsDegeneratedZone(aP2d, aS1, 2);
-    if (bFlag) {
-      return bFlag;
+
+    if(iter > 11) {
+      break;
     }
+    iter++;
   }
-  if (fabs(V1-VS12) < aDelta) {
-    bFlag=IsDegeneratedZone(aP2d, aS1, 2);
-    if (bFlag) {
-      return bFlag;
-    }
-  }
-  // Check on Surf 2
-  aD=aGAS2.UResolution(aDelta);
-  aP2d.SetCoord(U2, V2);
-  if (fabs(U2-US21) < aDelta) {
-    bFlag=IsDegeneratedZone(aP2d, aS2, 1);
-    if (bFlag) {
-      return bFlag;
-    }
-  }
-  if (fabs(U2-US22) < aDelta) {
-    bFlag=IsDegeneratedZone(aP2d, aS2, 1);
-    if (bFlag) {
-      return bFlag;
-    }
-  }
-  aD=aGAS2.VResolution(aDelta);
-  if (fabs(V2-VS21) < aDelta) {
-    bFlag=IsDegeneratedZone(aP2d, aS2, 2);
-    if (bFlag) {  
-      return bFlag;
-    }
-  }
-  if (fabs(V2-VS22) < aDelta) {
-    bFlag=IsDegeneratedZone(aP2d, aS2, 2);
-    if (bFlag) {
-      return bFlag;
-    }
-  }
-  return !bFlag;
-}
-
-//=======================================================================
-//function : IsDegeneratedZone
-//purpose  : 
-//=======================================================================
-Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
-                                   const Handle(Geom_Surface)& aS,
-                                   const Standard_Integer iDir)
-{
-  Standard_Boolean bFlag=Standard_True;
-  Standard_Real US1, US2, VS1, VS2, dY, dX, d1, d2, dD;
-  Standard_Real aXm, aYm, aXb, aYb, aXe, aYe;
-  aS->Bounds(US1, US2, VS1, VS2); 
-
-  gp_Pnt aPm, aPb, aPe;
-  
-  aXm=aP2d.X();
-  aYm=aP2d.Y();
-  
-  aS->D0(aXm, aYm, aPm); 
-  
-  dX=1.e-5;
-  dY=1.e-5;
-  dD=1.e-12;
-
-  if (iDir==1) {
-    aXb=aXm;
-    aXe=aXm;
-    aYb=aYm-dY;
-    if (aYb < VS1) {
-      aYb=VS1;
-    }
-    aYe=aYm+dY;
-    if (aYe > VS2) {
-      aYe=VS2;
-    }
-    aS->D0(aXb, aYb, aPb);
-    aS->D0(aXe, aYe, aPe);
-    
-    d1=aPm.Distance(aPb);
-    d2=aPm.Distance(aPe);
-    if (d1 < dD && d2 < dD) {
-      return bFlag;
-    }
-    return !bFlag;
-  }
-  //
-  else if (iDir==2) {
-    aYb=aYm;
-    aYe=aYm;
-    aXb=aXm-dX;
-    if (aXb < US1) {
-      aXb=US1;
-    }
-    aXe=aXm+dX;
-    if (aXe > US2) {
-      aXe=US2;
-    }
-    aS->D0(aXb, aYb, aPb);
-    aS->D0(aXe, aYe, aPe);
-    
-    d1=aPm.Distance(aPb);
-    d2=aPm.Distance(aPe);
-    if (d1 < dD && d2 < dD) {
-      return bFlag;
-    }
-    return !bFlag;
-  }
-  return !bFlag;
-}
-
-//=========================================================================
-// static function : ComputePurgedWLine
-// purpose : Removes equal points (leave one of equal points) from theWLine
-//           and recompute vertex parameters.
-//           Returns new WLine or null WLine if the number
-//           of the points is less than 2.
-//=========================================================================
-Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine) {
-  Standard_Integer i, k, v, nb, nbvtx;
-  Handle(IntPatch_WLine) aResult;
-  nbvtx = theWLine->NbVertex();
-  nb = theWLine->NbPnts();
-  if (nb==2) {
-    const IntSurf_PntOn2S& p1 = theWLine->Point(1);
-    const IntSurf_PntOn2S& p2 = theWLine->Point(2);
-    if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
-      return aResult;
-    }
-  }
-  //
-  Handle(IntPatch_WLine) aLocalWLine;
-  Handle(IntPatch_WLine) aTmpWLine = theWLine;
-  Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
-  aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
-  for(i = 1; i <= nb; i++) {
-    aLineOn2S->Add(theWLine->Point(i));
-  }
-
-  for(v = 1; v <= nbvtx; v++) {
-    aLocalWLine->AddVertex(theWLine->Vertex(v));
-  }
-  
-  for(i = 1; i <= aLineOn2S->NbPoints(); i++) {
-    Standard_Integer aStartIndex = i + 1;
-    Standard_Integer anEndIndex = i + 5;
-    nb = aLineOn2S->NbPoints();
-    anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
-
-    if((aStartIndex > nb) || (anEndIndex <= 1)) {
-      continue;
-    }
-    k = aStartIndex;
-
-    while(k <= anEndIndex) {
-      
-      if(i != k) {
-        IntSurf_PntOn2S p1 = aLineOn2S->Value(i);
-        IntSurf_PntOn2S p2 = aLineOn2S->Value(k);
-        
-        if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
-          aTmpWLine = aLocalWLine;
-          aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
-
-          for(v = 1; v <= aTmpWLine->NbVertex(); v++) {
-            IntPatch_Point aVertex = aTmpWLine->Vertex(v);
-            Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
-
-            if(avertexindex >= k) {
-              aVertex.SetParameter(aVertex.ParameterOnLine() - 1.);
-            }
-            aLocalWLine->AddVertex(aVertex);
-          }
-          aLineOn2S->RemovePoint(k);
-          anEndIndex--;
-          continue;
-        }
-      }
-      k++;
-    }
-  }
-
-  if(aLineOn2S->NbPoints() > 1) {
-    aResult = aLocalWLine;
-  }
-  return aResult;
-}
-
-//=======================================================================
-//function : TolR3d
-//purpose  : 
-//=======================================================================
-void TolR3d(const TopoDS_Face& aF1,
-            const TopoDS_Face& aF2,
-            Standard_Real& myTolReached3d)
-{
-  Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
-      
-  aTolTresh=2.999999e-3;
-  aTolF1 = BRep_Tool::Tolerance(aF1);
-  aTolF2 = BRep_Tool::Tolerance(aF2);
-  aTolFMax=Max(aTolF1, aTolF2);
-  
-  if (aTolFMax>aTolTresh) {
-    myTolReached3d=aTolFMax;
-  }
-}
-//=======================================================================
-//function : IsPointOnBoundary
-//purpose  : 
-//=======================================================================
-Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
-                                   const Standard_Real theFirstBoundary,
-                                   const Standard_Real theSecondBoundary,
-                                   const Standard_Real theResolution,
-                                   Standard_Boolean&   IsOnFirstBoundary) 
-{
-  Standard_Boolean bRet;
-  Standard_Integer i;
-  Standard_Real adist;
-  //
-  bRet=Standard_False;
-  for(i = 0; i < 2; ++i) {
-    IsOnFirstBoundary = (i == 0);
-    if (IsOnFirstBoundary) {
-      adist = fabs(theParameter - theFirstBoundary);
-    }
-    else {
-      adist = fabs(theParameter - theSecondBoundary);
-    }
-    if(adist < theResolution) {
-      return !bRet;
-    }
-  }
-  return bRet;
-}
-// ------------------------------------------------------------------------------------------------
-// static function: FindPoint
-// purpose:
-// ------------------------------------------------------------------------------------------------
-Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
-                           const gp_Pnt2d&     theLastPoint,
-                           const Standard_Real theUmin, 
-                           const Standard_Real theUmax,
-                           const Standard_Real theVmin,
-                           const Standard_Real theVmax,
-                           gp_Pnt2d&           theNewPoint) {
-  
-  gp_Vec2d aVec(theFirstPoint, theLastPoint);
-  Standard_Integer i = 0, j = 0;
-
-  for(i = 0; i < 4; i++) {
-    gp_Vec2d anOtherVec;
-    gp_Vec2d anOtherVecNormal;
-    gp_Pnt2d aprojpoint = theLastPoint;    
-
-    if((i % 2) == 0) {
-      anOtherVec.SetX(0.);
-      anOtherVec.SetY(1.);
-      anOtherVecNormal.SetX(1.);
-      anOtherVecNormal.SetY(0.);
-
-      if(i < 2)
-        aprojpoint.SetX(theUmin);
-      else
-        aprojpoint.SetX(theUmax);
-    }
-    else {
-      anOtherVec.SetX(1.);
-      anOtherVec.SetY(0.);
-      anOtherVecNormal.SetX(0.);
-      anOtherVecNormal.SetY(1.);
-
-      if(i < 2)
-        aprojpoint.SetY(theVmin);
-      else
-        aprojpoint.SetY(theVmax);
-    }
-    gp_Vec2d anormvec = aVec;
-    anormvec.Normalize();
-    RefineVector(anormvec);
-    Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
-
-    if(fabs(adot1) < Precision::Angular())
-      continue;
-    Standard_Real adist = 0.;
-    Standard_Boolean bIsOut = Standard_False;
-
-    if((i % 2) == 0) {
-      adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
-      bIsOut = (i < 2) ? (theLastPoint.X() < theUmin) : (theLastPoint.X() > theUmax);
-    }
-    else {
-      adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
-      bIsOut = (i < 2) ? (theLastPoint.Y() < theVmin) : (theLastPoint.Y() > theVmax);
-    }
-    Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
-
-    for(j = 0; j < 2; j++) {
-      anoffset = (j == 0) ? anoffset : -anoffset;
-      gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
-      gp_Vec2d acurvec(theLastPoint, acurpoint);
-      if ( bIsOut )
-        acurvec.Reverse();
-
-      Standard_Real aDotX, anAngleX;
-      //
-      aDotX = aVec.Dot(acurvec);
-      anAngleX = aVec.Angle(acurvec);
-      //
-      if(aDotX > 0. && fabs(anAngleX) < Precision::PConfusion()) {
-        if((i % 2) == 0) {
-          if((acurpoint.Y() >= theVmin) &&
-             (acurpoint.Y() <= theVmax)) {
-            theNewPoint = acurpoint;
-            return Standard_True;
-          }
-        }
-        else {
-          if((acurpoint.X() >= theUmin) &&
-             (acurpoint.X() <= theUmax)) {
-            theNewPoint = acurpoint;
-            return Standard_True;
-          }
-        }
-      }
-    }
-  }
-  return Standard_False;
-}
-
-
-// ------------------------------------------------------------------------------------------------
-// static function: FindPoint
-// purpose: Find point on the boundary of radial tangent zone
-// ------------------------------------------------------------------------------------------------
-Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
-                           const gp_Pnt2d&     theLastPoint,
-                           const Standard_Real theUmin, 
-                           const Standard_Real theUmax,
-                           const Standard_Real theVmin,
-                           const Standard_Real theVmax,
-                           const gp_Pnt2d&     theTanZoneCenter,
-                           const Standard_Real theZoneRadius,
-                           Handle(GeomAdaptor_HSurface) theGASurface,
-                           gp_Pnt2d&           theNewPoint) {
-  theNewPoint = theLastPoint;
-
-  if ( !IsInsideTanZone( theLastPoint, theTanZoneCenter, theZoneRadius, theGASurface) )
-    return Standard_False;
-
-  Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
-  Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
-
-  Standard_Real aRadius = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
-  gp_Ax22d anAxis( theTanZoneCenter, gp_Dir2d(1, 0), gp_Dir2d(0, 1) );
-  gp_Circ2d aCircle( anAxis, aRadius );
-  
-  //
-  gp_Vec2d aDir( theLastPoint.XY() - theFirstPoint.XY() );
-  Standard_Real aLength = aDir.Magnitude();
-  if ( aLength <= gp::Resolution() )
-    return Standard_False;
-  gp_Lin2d aLine( theFirstPoint, aDir );
-
-  //
-  Handle(Geom2d_Line) aCLine = new Geom2d_Line( aLine );
-  Handle(Geom2d_TrimmedCurve) aC1 = new Geom2d_TrimmedCurve( aCLine, 0, aLength );
-  Handle(Geom2d_Circle) aC2 = new Geom2d_Circle( aCircle );
-
-  Standard_Real aTol = aRadius * 0.001;
-  aTol = ( aTol < Precision::PConfusion() ) ? Precision::PConfusion() : aTol;
-
-  Geom2dAPI_InterCurveCurve anIntersector;
-  anIntersector.Init( aC1, aC2, aTol );
-
-  if ( anIntersector.NbPoints() == 0 )
-    return Standard_False;
-
-  Standard_Boolean aFound = Standard_False;
-  Standard_Real aMinDist = aLength * aLength;
-  Standard_Integer i = 0;
-  for ( i = 1; i <= anIntersector.NbPoints(); i++ ) {
-    gp_Pnt2d aPInt = anIntersector.Point( i );
-    if ( aPInt.SquareDistance( theFirstPoint ) < aMinDist ) {
-      if ( ( aPInt.X() >= theUmin ) && ( aPInt.X() <= theUmax ) &&
-           ( aPInt.Y() >= theVmin ) && ( aPInt.Y() <= theVmax ) ) {
-        theNewPoint = aPInt;
-        aFound = Standard_True;
-      }
-    }
-  }
-
-  return aFound;
-}
-
-// ------------------------------------------------------------------------------------------------
-// static function: IsInsideTanZone
-// purpose: Check if point is inside a radial tangent zone
-// ------------------------------------------------------------------------------------------------
-Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
-                                 const gp_Pnt2d&     theTanZoneCenter,
-                                 const Standard_Real theZoneRadius,
-                                 Handle(GeomAdaptor_HSurface) theGASurface) {
-
-  Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
-  Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
-  Standard_Real aRadiusSQR = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
-  aRadiusSQR *= aRadiusSQR;
-  if ( thePoint.SquareDistance( theTanZoneCenter ) <= aRadiusSQR )
-    return Standard_True;
-  return Standard_False;
-}
-
-// ------------------------------------------------------------------------------------------------
-// static function: CheckTangentZonesExist
-// purpose: Check if tangent zone exists
-// ------------------------------------------------------------------------------------------------
-Standard_Boolean CheckTangentZonesExist( const Handle(GeomAdaptor_HSurface)& theSurface1,
-                                        const Handle(GeomAdaptor_HSurface)&  theSurface2 ) 
-{
-  if ( ( theSurface1->GetType() != GeomAbs_Torus ) ||
-      ( theSurface2->GetType() != GeomAbs_Torus ) )
-    return Standard_False;
-
-  gp_Torus aTor1 = theSurface1->Torus();
-  gp_Torus aTor2 = theSurface2->Torus();
-
-  if ( aTor1.Location().Distance( aTor2.Location() ) > Precision::Confusion() )
-    return Standard_False;
-
-  if ( ( fabs( aTor1.MajorRadius() - aTor2.MajorRadius() ) > Precision::Confusion() ) ||
-       ( fabs( aTor1.MinorRadius() - aTor2.MinorRadius() ) > Precision::Confusion() ) )
-    return Standard_False;
-
-  if ( ( aTor1.MajorRadius() < aTor1.MinorRadius() ) ||
-       ( aTor2.MajorRadius() < aTor2.MinorRadius() ) )
-    return Standard_False;
-  return Standard_True;
-}
-
-// ------------------------------------------------------------------------------------------------
-// static function: ComputeTangentZones
-// purpose: 
-// ------------------------------------------------------------------------------------------------
-Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
-                                     const Handle(GeomAdaptor_HSurface)&  theSurface2,
-                                     const TopoDS_Face&                   theFace1,
-                                     const TopoDS_Face&                   theFace2,
-                                     Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
-                                     Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
-                                     Handle(TColStd_HArray1OfReal)&       theResultRadius,
-                                     const Handle(IntTools_Context)& aContext)
-{
-  Standard_Integer aResult = 0;
-  if ( !CheckTangentZonesExist( theSurface1, theSurface2 ) )
-    return aResult;
-
-
-  TColgp_SequenceOfPnt2d aSeqResultS1, aSeqResultS2;
-  TColStd_SequenceOfReal aSeqResultRad;
-
-  gp_Torus aTor1 = theSurface1->Torus();
-  gp_Torus aTor2 = theSurface2->Torus();
-
-  gp_Ax2 anax1( aTor1.Location(), aTor1.Axis().Direction() );
-  gp_Ax2 anax2( aTor2.Location(), aTor2.Axis().Direction() );
-  Standard_Integer j = 0;
-
-  for ( j = 0; j < 2; j++ ) {
-    Standard_Real aCoef = ( j == 0 ) ? -1 : 1;
-    Standard_Real aRadius1 = fabs(aTor1.MajorRadius() + aCoef * aTor1.MinorRadius());
-    Standard_Real aRadius2 = fabs(aTor2.MajorRadius() + aCoef * aTor2.MinorRadius());
-
-    gp_Circ aCircle1( anax1, aRadius1 );
-    gp_Circ aCircle2( anax2, aRadius2 );
-
-    // roughly compute radius of tangent zone for perpendicular case
-    Standard_Real aCriteria = Precision::Confusion() * 0.5;
-
-    Standard_Real aT1 = aCriteria;
-    Standard_Real aT2 = aCriteria;
-    if ( j == 0 ) {
-      // internal tangency
-      Standard_Real aR = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
-      //aT1 = aCriteria * aCriteria + aR * aR - ( aR - aCriteria ) * ( aR - aCriteria );
-      aT1 = 2. * aR * aCriteria;
-      aT2 = aT1;
-    }
-    else {
-      // external tangency
-      Standard_Real aRb = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
-      Standard_Real aRm = ( aRadius1 < aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
-      Standard_Real aDelta = aRb - aCriteria;
-      aDelta *= aDelta;
-      aDelta -= aRm * aRm;
-      aDelta /= 2. * (aRb - aRm);
-      aDelta -= 0.5 * (aRb - aRm);
-      
-      aT1 = 2. * aRm * (aRm - aDelta);
-      aT2 = aT1;
-    }
-    aCriteria = ( aT1 > aT2) ? aT1 : aT2;
-    if ( aCriteria > 0 )
-      aCriteria = sqrt( aCriteria );
-
-    if ( aCriteria > 0.5 * aTor1.MinorRadius() ) {
-      // too big zone -> drop to minimum
-      aCriteria = Precision::Confusion();
-    }
-
-    GeomAdaptor_Curve aC1( new Geom_Circle(aCircle1) );
-    GeomAdaptor_Curve aC2( new Geom_Circle(aCircle2) );
-    Extrema_ExtCC anExtrema(aC1, aC2, 0, 2. * M_PI, 0, 2. * M_PI, 
-                            Precision::PConfusion(), Precision::PConfusion());
-            
-    if ( anExtrema.IsDone() ) {
-
-      Standard_Integer i = 0;
-      for ( i = 1; i <= anExtrema.NbExt(); i++ ) {
-        if ( anExtrema.SquareDistance(i) > aCriteria * aCriteria )
-          continue;
-
-        Extrema_POnCurv P1, P2;
-        anExtrema.Points( i, P1, P2 );
-
-        Standard_Boolean bFoundResult = Standard_True;
-        gp_Pnt2d pr1, pr2;
-
-        Standard_Integer surfit = 0;
-        for ( surfit = 0; surfit < 2; surfit++ ) {
-          GeomAPI_ProjectPointOnSurf& aProjector = 
-            (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
-
-          gp_Pnt aP3d = (surfit == 0) ? P1.Value() : P2.Value();
-          aProjector.Perform(aP3d);
-
-          if(!aProjector.IsDone())
-            bFoundResult = Standard_False;
-          else {
-            if(aProjector.LowerDistance() > aCriteria) {
-              bFoundResult = Standard_False;
-            }
-            else {
-              Standard_Real foundU = 0, foundV = 0;
-              aProjector.LowerDistanceParameters(foundU, foundV);
-              if ( surfit == 0 )
-                pr1 = gp_Pnt2d( foundU, foundV );
-              else
-                pr2 = gp_Pnt2d( foundU, foundV );
-            }
-          }
-        }
-        if ( bFoundResult ) {
-          aSeqResultS1.Append( pr1 );
-          aSeqResultS2.Append( pr2 );
-          aSeqResultRad.Append( aCriteria );
-
-          // torus is u and v periodic
-          const Standard_Real twoPI = M_PI + M_PI;
-          Standard_Real arr1tmp[2] = {pr1.X(), pr1.Y()};
-          Standard_Real arr2tmp[2] = {pr2.X(), pr2.Y()};
-
-          // iteration on period bounds
-          for ( Standard_Integer k1 = 0; k1 < 2; k1++ ) {
-            Standard_Real aBound = ( k1 == 0 ) ? 0 : twoPI;
-            Standard_Real aShift = ( k1 == 0 ) ? twoPI : -twoPI;
-
-            // iteration on surfaces
-            for ( Standard_Integer k2 = 0; k2 < 2; k2++ ) {
-              Standard_Real* arr1 = ( k2 == 0 ) ? arr1tmp : arr2tmp;
-              Standard_Real* arr2 = ( k2 != 0 ) ? arr1tmp : arr2tmp;
-              TColgp_SequenceOfPnt2d& aSeqS1 = ( k2 == 0 ) ? aSeqResultS1 : aSeqResultS2; 
-              TColgp_SequenceOfPnt2d& aSeqS2 = ( k2 != 0 ) ? aSeqResultS1 : aSeqResultS2; 
-
-              if (fabs(arr1[0] - aBound) < Precision::PConfusion()) {
-                aSeqS1.Append( gp_Pnt2d( arr1[0] + aShift, arr1[1] ) );
-                aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
-                aSeqResultRad.Append( aCriteria );
-              }
-              if (fabs(arr1[1] - aBound) < Precision::PConfusion()) {
-                aSeqS1.Append( gp_Pnt2d( arr1[0], arr1[1] + aShift) );
-                aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
-                aSeqResultRad.Append( aCriteria );
-              }
-            }
-          } //
-        }
-      }
-    }
-  }
-  aResult = aSeqResultRad.Length();
-
-  if ( aResult > 0 ) {
-    theResultOnS1 = new TColgp_HArray1OfPnt2d( 1, aResult );
-    theResultOnS2 = new TColgp_HArray1OfPnt2d( 1, aResult );
-    theResultRadius = new TColStd_HArray1OfReal( 1, aResult );
-
-    for ( Standard_Integer i = 1 ; i <= aResult; i++ ) {
-      theResultOnS1->SetValue( i, aSeqResultS1.Value(i) );
-      theResultOnS2->SetValue( i, aSeqResultS2.Value(i) );
-      theResultRadius->SetValue( i, aSeqResultRad.Value(i) );
-    }
-  }
-  return aResult;
-}
-
-// ------------------------------------------------------------------------------------------------
-// static function: AdjustByNeighbour
-// purpose:
-// ------------------------------------------------------------------------------------------------
-gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
-                           const gp_Pnt2d&     theOriginalPoint,
-                           Handle(GeomAdaptor_HSurface) theGASurface) {
-  
-  gp_Pnt2d ap1 = theaNeighbourPoint;
-  gp_Pnt2d ap2 = theOriginalPoint;
-
-  if ( theGASurface->IsUPeriodic() ) {
-    Standard_Real aPeriod     = theGASurface->UPeriod();
-    gp_Pnt2d aPTest = ap2;
-    Standard_Real aSqDistMin = 1.e+100;
-
-    for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
-      aPTest.SetX( theOriginalPoint.X() + aPeriod * pIt );
-      Standard_Real dd = ap1.SquareDistance( aPTest );
-
-      if ( dd < aSqDistMin ) {
-        ap2 = aPTest;
-        aSqDistMin = dd;
-      }
-    }
-  }
-  if ( theGASurface->IsVPeriodic() ) {
-    Standard_Real aPeriod     = theGASurface->VPeriod();
-    gp_Pnt2d aPTest = ap2;
-    Standard_Real aSqDistMin = 1.e+100;
-
-    for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
-      aPTest.SetY( theOriginalPoint.Y() + aPeriod * pIt );
-      Standard_Real dd = ap1.SquareDistance( aPTest );
-
-      if ( dd < aSqDistMin ) {
-        ap2 = aPTest;
-        aSqDistMin = dd;
-      }
-    }
-  }
-  return ap2;
-}
-
-// ------------------------------------------------------------------------------------------------
-//function: DecompositionOfWLine
-// purpose:
-// ------------------------------------------------------------------------------------------------
-Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
-                                      const Handle(GeomAdaptor_HSurface)&            theSurface1, 
-                                      const Handle(GeomAdaptor_HSurface)&            theSurface2,
-                                      const TopoDS_Face&                             theFace1,
-                                      const TopoDS_Face&                             theFace2,
-                                      const GeomInt_LineConstructor&                 theLConstructor,
-                                      const Standard_Boolean                         theAvoidLConstructor,
-                                      IntPatch_SequenceOfLine&                       theNewLines,
-                                      Standard_Real&                                 theReachedTol3d,
-                                      const Handle(IntTools_Context)& aContext) 
-{
-
-  Standard_Boolean bRet, bAvoidLineConstructor;
-  Standard_Integer aNbPnts, aNbParts;
-  //
-  bRet=Standard_False;
-  aNbPnts=theWLine->NbPnts();
-  bAvoidLineConstructor=theAvoidLConstructor;
-  //
-  if(!aNbPnts){
-    return bRet;
-  }
-  if (!bAvoidLineConstructor) {
-    aNbParts=theLConstructor.NbParts();
-    if (!aNbParts) {
-      return bRet;
-    }
-  }
-  //
-  Standard_Boolean bIsPrevPointOnBoundary, bIsPointOnBoundary, bIsCurrentPointOnBoundary;
-  Standard_Integer nblines, pit, i, j;
-  Standard_Real aTol;
-  TColStd_Array1OfListOfInteger anArrayOfLines(1, aNbPnts); 
-  TColStd_Array1OfInteger       anArrayOfLineType(1, aNbPnts);
-  TColStd_ListOfInteger aListOfPointIndex;
-  
-  Handle(TColgp_HArray1OfPnt2d) aTanZoneS1;
-  Handle(TColgp_HArray1OfPnt2d) aTanZoneS2;
-  Handle(TColStd_HArray1OfReal) aTanZoneRadius;
-  Standard_Integer aNbZone = ComputeTangentZones( theSurface1, theSurface2, theFace1, theFace2,
-                                                 aTanZoneS1, aTanZoneS2, aTanZoneRadius, aContext);
-  
-  //
-  nblines=0;
-  aTol=Precision::Confusion();
-  aTol=0.5*aTol;
-  bIsPrevPointOnBoundary=Standard_False;
-  bIsPointOnBoundary=Standard_False;
-  //
-  // 1. ...
-  //
-  // Points
-  for(pit = 1; pit <= aNbPnts; ++pit) {
-    Standard_Boolean bIsOnFirstBoundary, isperiodic;
-    Standard_Real aResolution, aPeriod, alowerboundary, aupperboundary, U, V;
-    Standard_Real aParameter, anoffset, anAdjustPar;
-    Standard_Real umin, umax, vmin, vmax;
-    //
-    bIsCurrentPointOnBoundary = Standard_False;
-    const IntSurf_PntOn2S& aPoint = theWLine->Point(pit);
-    //
-    // Surface
-    for(i = 0; i < 2; ++i) {
-      Handle(GeomAdaptor_HSurface) aGASurface = (!i) ? theSurface1 : theSurface2;
-      aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
-      if(!i) {
-        aPoint.ParametersOnS1(U, V);
-      }
-      else {
-        aPoint.ParametersOnS2(U, V);
-      }
-      // U, V
-      for(j = 0; j < 2; j++) {
-        isperiodic = (!j) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
-        if(!isperiodic){
-          continue;
-        }
-        //
-        if (!j) {
-          aResolution=aGASurface->UResolution(aTol);
-          aPeriod=aGASurface->UPeriod();
-          alowerboundary=umin;
-          aupperboundary=umax;
-          aParameter=U;
-        }
-        else {
-          aResolution=aGASurface->VResolution(aTol);
-          aPeriod=aGASurface->VPeriod();
-          alowerboundary=vmin;
-          aupperboundary=vmax;
-          aParameter=V;
-        }
-        
-        GeomInt::AdjustPeriodic(aParameter, 
-                                       alowerboundary, 
-                                       aupperboundary, 
-                                       aPeriod,
-                                       anAdjustPar,
-                                       anoffset);
-        //
-        bIsOnFirstBoundary = Standard_True;// ?
-        bIsPointOnBoundary=
-          IsPointOnBoundary(anAdjustPar, 
-                            alowerboundary, 
-                            aupperboundary,
-                            aResolution, 
-                            bIsOnFirstBoundary);
-        //
-        if(bIsPointOnBoundary) {
-          bIsCurrentPointOnBoundary = Standard_True;
-          break;
-        }
-        else {
-          // check if a point belong to a tangent zone. Begin
-          Standard_Integer zIt = 0;
-          for ( zIt = 1; zIt <= aNbZone; zIt++ ) {
-            gp_Pnt2d aPZone = (i == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
-            Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
-
-            if ( IsInsideTanZone(gp_Pnt2d( U, V ), aPZone, aZoneRadius, aGASurface ) ) {
-              // set boundary flag to split the curve by a tangent zone
-              bIsPointOnBoundary = Standard_True;
-              bIsCurrentPointOnBoundary = Standard_True;
-              if ( theReachedTol3d < aZoneRadius ) {
-                theReachedTol3d = aZoneRadius;
-              }
-              break;
-            }
-          }
-        }
-      }//for(j = 0; j < 2; j++) {
-
-      if(bIsCurrentPointOnBoundary){
-        break;
-      }
-    }//for(i = 0; i < 2; ++i) {
-    //
-    if((bIsCurrentPointOnBoundary != bIsPrevPointOnBoundary)) {
-      if(!aListOfPointIndex.IsEmpty()) {
-        nblines++;
-        anArrayOfLines.SetValue(nblines, aListOfPointIndex);
-        anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
-        aListOfPointIndex.Clear();
-      }
-      bIsPrevPointOnBoundary = bIsCurrentPointOnBoundary;
-    }
-    aListOfPointIndex.Append(pit);
-  } //for(pit = 1; pit <= aNbPnts; ++pit) {
-  //
-  if(!aListOfPointIndex.IsEmpty()) {
-    nblines++;
-    anArrayOfLines.SetValue(nblines, aListOfPointIndex);
-    anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
-    aListOfPointIndex.Clear();
-  }
-  //
-  if(nblines<=1) {
-    return bRet; //Standard_False;
-  }
-  //
-  // 
-  // 2. Correct wlines.begin
-  TColStd_Array1OfListOfInteger anArrayOfLineEnds(1, nblines);
-  Handle(IntSurf_LineOn2S) aSeqOfPntOn2S = new IntSurf_LineOn2S();
-  //
-  for(i = 1; i <= nblines; i++) {
-    if(anArrayOfLineType.Value(i) != 0) {
-      continue;
-    }
-    const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
-    TColStd_ListOfInteger aListOfFLIndex;
-
-    for(j = 0; j < 2; j++) {
-      Standard_Integer aneighbourindex = (j == 0) ? (i - 1) : (i + 1);
-
-      if((aneighbourindex < 1) || (aneighbourindex > nblines))
-        continue;
-
-      if(anArrayOfLineType.Value(aneighbourindex) == 0)
-        continue;
-      const TColStd_ListOfInteger& aNeighbour = anArrayOfLines.Value(aneighbourindex);
-      Standard_Integer anIndex = (j == 0) ? aNeighbour.Last() : aNeighbour.First();
-      const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
-
-      IntSurf_PntOn2S aNewP = aPoint;
-      if(aListOfIndex.Extent() < 2) {
-        aSeqOfPntOn2S->Add(aNewP);
-        aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
-        continue;
-      }
-      //
-      Standard_Integer iFirst = aListOfIndex.First();
-      Standard_Integer iLast  = aListOfIndex.Last();
-      //
-      for(Standard_Integer surfit = 0; surfit < 2; surfit++) {
-
-        Handle(GeomAdaptor_HSurface) aGASurface = (surfit == 0) ? theSurface1 : theSurface2;
-        Standard_Real umin=0., umax=0., vmin=0., vmax=0.;
-        aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
-        Standard_Real U=0., V=0.;
-
-        if(surfit == 0)
-          aNewP.ParametersOnS1(U, V);
-        else
-          aNewP.ParametersOnS2(U, V);
-        Standard_Integer nbboundaries = 0;
-
-        Standard_Boolean bIsNearBoundary = Standard_False;
-        Standard_Integer aZoneIndex = 0;
-        Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
-        Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
-        
-
-        for(Standard_Integer parit = 0; parit < 2; parit++) {
-          Standard_Boolean isperiodic = (parit == 0) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
-
-          Standard_Real aResolution = (parit == 0) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
-          Standard_Real alowerboundary = (parit == 0) ? umin : vmin;
-          Standard_Real aupperboundary = (parit == 0) ? umax : vmax;
-
-          Standard_Real aParameter = (parit == 0) ? U : V;
-          Standard_Boolean bIsOnFirstBoundary = Standard_True;
-  
-          if(!isperiodic) {
-            bIsPointOnBoundary=
-              IsPointOnBoundary(aParameter, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
-            if(bIsPointOnBoundary) {
-              bIsUBoundary = (parit == 0);
-              bIsFirstBoundary = bIsOnFirstBoundary;
-              nbboundaries++;
-            }
-          }
-          else {
-            Standard_Real aPeriod     = (parit == 0) ? aGASurface->UPeriod() : aGASurface->VPeriod();
-            Standard_Real anoffset, anAdjustPar;
-            GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
-                                           aPeriod, anAdjustPar, anoffset);
-
-            bIsPointOnBoundary=
-              IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
-            if(bIsPointOnBoundary) {
-              bIsUBoundary = (parit == 0);
-              bIsFirstBoundary = bIsOnFirstBoundary;
-              nbboundaries++;
-            }
-            else {
-            //check neighbourhood of boundary
-            Standard_Real anEpsilon = aResolution * 100.;
-            Standard_Real aPart = ( aupperboundary - alowerboundary ) * 0.1;
-            anEpsilon = ( anEpsilon > aPart ) ? aPart : anEpsilon;
-            
-              bIsNearBoundary = IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, 
-                                                  anEpsilon, bIsOnFirstBoundary);
-
-            }
-          }
-        }
-
-        // check if a point belong to a tangent zone. Begin
-        for ( Standard_Integer zIt = 1; zIt <= aNbZone; zIt++ ) {
-          gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
-          Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
-
-          Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
-          const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
-          Standard_Real nU1, nV1;
-            
-          if(surfit == 0)
-            aNeighbourPoint.ParametersOnS1(nU1, nV1);
-          else
-            aNeighbourPoint.ParametersOnS2(nU1, nV1);
-          gp_Pnt2d ap1(nU1, nV1);
-          gp_Pnt2d ap2 = AdjustByNeighbour( ap1, gp_Pnt2d( U, V ), aGASurface );
-
-
-          if ( IsInsideTanZone( ap2, aPZone, aZoneRadius, aGASurface ) ) {
-            aZoneIndex = zIt;
-            bIsNearBoundary = Standard_True;
-            if ( theReachedTol3d < aZoneRadius ) {
-              theReachedTol3d = aZoneRadius;
-            }
-          }
-        }
-        // check if a point belong to a tangent zone. End
-        Standard_Boolean bComputeLineEnd = Standard_False;
-
-        if(nbboundaries == 2) {
-          //xf
-          bComputeLineEnd = Standard_True;
-          //xt
-        }
-        else if(nbboundaries == 1) {
-          Standard_Boolean isperiodic = (bIsUBoundary) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
-
-          if(isperiodic) {
-            Standard_Real alowerboundary = (bIsUBoundary) ? umin : vmin;
-            Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
-            Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
-            Standard_Real aParameter = (bIsUBoundary) ? U : V;
-            Standard_Real anoffset, anAdjustPar;
-            GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary, 
-                                           aPeriod, anAdjustPar, anoffset);
-
-            Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
-            Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
-            anotherPar += anoffset;
-            Standard_Integer aneighbourpointindex = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
-            const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex);
-            Standard_Real nU1, nV1;
-
-            if(surfit == 0)
-              aNeighbourPoint.ParametersOnS1(nU1, nV1);
-            else
-              aNeighbourPoint.ParametersOnS2(nU1, nV1);
-            
-            Standard_Real adist1 = (bIsUBoundary) ? fabs(nU1 - U) : fabs(nV1 - V);
-            Standard_Real adist2 = (bIsUBoundary) ? fabs(nU1 - anotherPar) : fabs(nV1 - anotherPar);
-            bComputeLineEnd = Standard_True;
-            Standard_Boolean bCheckAngle1 = Standard_False;
-            Standard_Boolean bCheckAngle2 = Standard_False;
-            gp_Vec2d aNewVec;
-            Standard_Real anewU = (bIsUBoundary) ? anotherPar : U;
-            Standard_Real anewV = (bIsUBoundary) ? V : anotherPar;
-
-            if(((adist1 - adist2) > Precision::PConfusion()) && 
-               (adist2 < (aPeriod / 4.))) {
-              bCheckAngle1 = Standard_True;
-              aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
-
-              if(aNewVec.SquareMagnitude() < gp::Resolution()) {
-                aNewP.SetValue((surfit == 0), anewU, anewV);
-                bCheckAngle1 = Standard_False;
-              }
-            }
-            else if(adist1 < (aPeriod / 4.)) {
-              bCheckAngle2 = Standard_True;
-              aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
-
-              if(aNewVec.SquareMagnitude() < gp::Resolution()) {
-                bCheckAngle2 = Standard_False;
-              }
-            }
-
-            if(bCheckAngle1 || bCheckAngle2) {
-              // assume there are at least two points in line (see "if" above)
-              Standard_Integer anindexother = aneighbourpointindex;
-
-              while((anindexother <= iLast) && (anindexother >= iFirst)) {
-                anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
-                const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
-                Standard_Real nU2, nV2;
-                
-                if(surfit == 0)
-                  aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
-                else
-                  aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
-                gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
-
-                if(aVecOld.SquareMagnitude() <= gp::Resolution()) {
-                  continue;
-                }
-                else {
-                  Standard_Real anAngle = aNewVec.Angle(aVecOld);
-
-                  if((fabs(anAngle) < (M_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
-
-                    if(bCheckAngle1) {
-                      Standard_Real U1, U2, V1, V2;
-                      IntSurf_PntOn2S atmppoint = aNewP;
-                      atmppoint.SetValue((surfit == 0), anewU, anewV);
-                      atmppoint.Parameters(U1, V1, U2, V2);
-                      gp_Pnt P1 = theSurface1->Value(U1, V1);
-                      gp_Pnt P2 = theSurface2->Value(U2, V2);
-                      gp_Pnt P0 = aPoint.Value();
-
-                      if(P0.IsEqual(P1, aTol) &&
-                         P0.IsEqual(P2, aTol) &&
-                         P1.IsEqual(P2, aTol)) {
-                        bComputeLineEnd = Standard_False;
-                        aNewP.SetValue((surfit == 0), anewU, anewV);
-                      }
-                    }
-
-                    if(bCheckAngle2) {
-                      bComputeLineEnd = Standard_False;
-                    }
-                  }
-                  break;
-                }
-              } // end while(anindexother...)
-            }
-          }
-        }
-        else if ( bIsNearBoundary ) {
-          bComputeLineEnd = Standard_True;
-        }
-
-        if(bComputeLineEnd) {
-
-          gp_Pnt2d anewpoint;
-          Standard_Boolean found = Standard_False;
-
-          if ( bIsNearBoundary ) {
-            // re-compute point near natural boundary or near tangent zone
-            Standard_Real u1, v1, u2, v2;
-            aNewP.Parameters( u1, v1, u2, v2 );
-            if(surfit == 0)
-              anewpoint = gp_Pnt2d( u1, v1 );
-            else
-              anewpoint = gp_Pnt2d( u2, v2 );
-            
-            Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
-            const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
-            Standard_Real nU1, nV1;
-            
-            if(surfit == 0)
-              aNeighbourPoint.ParametersOnS1(nU1, nV1);
-            else
-              aNeighbourPoint.ParametersOnS2(nU1, nV1);
-            gp_Pnt2d ap1(nU1, nV1);
-            gp_Pnt2d ap2;
-
-
-            if ( aZoneIndex ) {
-              // exclude point from a tangent zone
-              anewpoint = AdjustByNeighbour( ap1, anewpoint, aGASurface );
-              gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(aZoneIndex) : aTanZoneS2->Value(aZoneIndex);
-              Standard_Real aZoneRadius = aTanZoneRadius->Value(aZoneIndex);
-
-              if ( FindPoint(ap1, anewpoint, umin, umax, vmin, vmax, 
-                             aPZone, aZoneRadius, aGASurface, ap2) ) {
-                anewpoint = ap2;
-                found = Standard_True;
-              }
-            }
-            else if ( aGASurface->IsUPeriodic() || aGASurface->IsVPeriodic() ) {
-              // re-compute point near boundary if shifted on a period
-              ap2 = AdjustByNeighbour( ap1, anewpoint, aGASurface );
-
-              if ( ( ap2.X() < umin ) || ( ap2.X() > umax ) ||
-                  ( ap2.Y() < vmin ) || ( ap2.Y() > vmax ) ) {
-                found = FindPoint(ap1, ap2, umin, umax, vmin, vmax, anewpoint);
-              }
-              else {
-                anewpoint = ap2;
-                aNewP.SetValue( (surfit == 0), anewpoint.X(), anewpoint.Y() );
-              }
-            }
-          }
-          else {
-
-            Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
-            const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
-            Standard_Real nU1, nV1;
-
-            if(surfit == 0)
-              aNeighbourPoint.ParametersOnS1(nU1, nV1);
-            else
-              aNeighbourPoint.ParametersOnS2(nU1, nV1);
-            gp_Pnt2d ap1(nU1, nV1);
-            gp_Pnt2d ap2(nU1, nV1);
-            Standard_Integer aneighbourpointindex2 = aneighbourpointindex1;
-
-            while((aneighbourpointindex2 <= iLast) && (aneighbourpointindex2 >= iFirst)) {
-              aneighbourpointindex2 = (j == 0) ? (aneighbourpointindex2 + 1) : (aneighbourpointindex2 - 1);
-              const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(aneighbourpointindex2);
-              Standard_Real nU2, nV2;
-
-              if(surfit == 0)
-                aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
-              else
-                aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
-              ap2.SetX(nU2);
-              ap2.SetY(nV2);
-
-              if(ap1.SquareDistance(ap2) > gp::Resolution()) {
-                break;
-              }
-            }  
-            found = FindPoint(ap2, ap1, umin, umax, vmin, vmax, anewpoint);
-          }
-
-          if(found) {
-            // check point
-            Standard_Real aCriteria = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
-            GeomAPI_ProjectPointOnSurf& aProjector = 
-              (surfit == 0) ? aContext->ProjPS(theFace2) : aContext->ProjPS(theFace1);
-            Handle(GeomAdaptor_HSurface) aSurface = (surfit == 0) ? theSurface1 : theSurface2;
-
-            Handle(GeomAdaptor_HSurface) aSurfaceOther = (surfit == 0) ? theSurface2 : theSurface1;
-
-            gp_Pnt aP3d = aSurface->Value(anewpoint.X(), anewpoint.Y());
-            aProjector.Perform(aP3d);
-
-            if(aProjector.IsDone()) {
-              if(aProjector.LowerDistance() < aCriteria) {
-                Standard_Real foundU = U, foundV = V;
-                aProjector.LowerDistanceParameters(foundU, foundV);
-
-                //Correction of projected coordinates. Begin
-                //Note, it may be shifted on a period
-                Standard_Integer aneindex1 = (j == 0) ? iFirst : iLast;
-                const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneindex1);
-                Standard_Real nUn, nVn;
-                
-                if(surfit == 0)
-                  aNeighbourPoint.ParametersOnS2(nUn, nVn);
-                else
-                  aNeighbourPoint.ParametersOnS1(nUn, nVn);
-                gp_Pnt2d aNeighbour2d(nUn, nVn);
-                gp_Pnt2d anAdjustedPoint = AdjustByNeighbour( aNeighbour2d, gp_Pnt2d(foundU, foundV), aSurfaceOther );
-                foundU = anAdjustedPoint.X();
-                foundV = anAdjustedPoint.Y();
-
-                if ( ( anAdjustedPoint.X() < umin ) && ( anAdjustedPoint.X() > umax ) &&
-                    ( anAdjustedPoint.Y() < vmin ) && ( anAdjustedPoint.Y() > vmax ) ) {
-                  // attempt to roughly re-compute point
-                  foundU = ( foundU < umin ) ? umin : foundU;
-                  foundU = ( foundU > umax ) ? umax : foundU;
-                  foundV = ( foundV < vmin ) ? vmin : foundV;
-                  foundV = ( foundV > vmax ) ? vmax : foundV;
-
-                  GeomAPI_ProjectPointOnSurf& aProjector2 = 
-                    (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
-
-                  aP3d = aSurfaceOther->Value(foundU, foundV);
-                  aProjector2.Perform(aP3d);
-                  
-                  if(aProjector2.IsDone()) {
-                    if(aProjector2.LowerDistance() < aCriteria) {
-                      Standard_Real foundU2 = anewpoint.X(), foundV2 = anewpoint.Y();
-                      aProjector2.LowerDistanceParameters(foundU2, foundV2);
-                      anewpoint.SetX(foundU2);
-                      anewpoint.SetY(foundV2);
-                    }
-                  }
-                }
-                //Correction of projected coordinates. End
-
-                if(surfit == 0)
-                  aNewP.SetValue(aP3d, anewpoint.X(), anewpoint.Y(), foundU, foundV);
-                else
-                  aNewP.SetValue(aP3d, foundU, foundV, anewpoint.X(), anewpoint.Y());
-              }
-            }
-          }
-        }
-      }
-      aSeqOfPntOn2S->Add(aNewP);
-      aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
-    }
-    anArrayOfLineEnds.SetValue(i, aListOfFLIndex);
-  }
-  // Correct wlines.end
-
-  // Split wlines.begin
-  Standard_Integer nbiter;
-  //
-  nbiter=1;
-  if (!bAvoidLineConstructor) {
-    nbiter=theLConstructor.NbParts();
-  }
-  //
-  for(j = 1; j <= nbiter; ++j) {
-    Standard_Real fprm, lprm;
-    Standard_Integer ifprm, ilprm;
-    //
-    if(bAvoidLineConstructor) {
-      ifprm = 1;
-      ilprm = theWLine->NbPnts();
-    }
-    else {
-      theLConstructor.Part(j, fprm, lprm);
-      ifprm = (Standard_Integer)fprm;
-      ilprm = (Standard_Integer)lprm;
-    }
-
-    Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
-    //
-    for(i = 1; i <= nblines; i++) {
-      if(anArrayOfLineType.Value(i) != 0) {
-        continue;
-      }
-      const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
-      const TColStd_ListOfInteger& aListOfFLIndex = anArrayOfLineEnds.Value(i);
-      Standard_Boolean bhasfirstpoint = (aListOfFLIndex.Extent() == 2);
-      Standard_Boolean bhaslastpoint = (aListOfFLIndex.Extent() == 2);
-
-      if(!bhasfirstpoint && !aListOfFLIndex.IsEmpty()) {
-        bhasfirstpoint = (i != 1);
-      }
-
-      if(!bhaslastpoint && !aListOfFLIndex.IsEmpty()) {
-        bhaslastpoint = (i != nblines);
-      }
-      
-      Standard_Integer iFirst = aListOfIndex.First();
-      Standard_Integer iLast  = aListOfIndex.Last();
-      Standard_Boolean bIsFirstInside = ((ifprm >= iFirst) && (ifprm <= iLast));
-      Standard_Boolean bIsLastInside =  ((ilprm >= iFirst) && (ilprm <= iLast));
-
-      if(!bIsFirstInside && !bIsLastInside) {
-        if((ifprm < iFirst) && (ilprm > iLast)) {
-          // append whole line, and boundaries if neccesary
-          if(bhasfirstpoint) {
-            pit = aListOfFLIndex.First();
-            const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
-            aLineOn2S->Add(aP);
-          }
-          TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
-
-          for(; anIt.More(); anIt.Next()) {
-            pit = anIt.Value();
-            const IntSurf_PntOn2S& aP = theWLine->Point(pit);
-            aLineOn2S->Add(aP);
-          }
-
-          if(bhaslastpoint) {
-            pit = aListOfFLIndex.Last();
-            const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
-            aLineOn2S->Add(aP);
-          }
-
-          // check end of split line (end is almost always)
-          Standard_Integer aneighbour = i + 1;
-          Standard_Boolean bIsEndOfLine = Standard_True;
-
-          if(aneighbour <= nblines) {
-            const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
-
-            if((anArrayOfLineType.Value(aneighbour) != 0) &&
-               (aListOfNeighbourIndex.IsEmpty())) {
-              bIsEndOfLine = Standard_False;
-            }
-          }
-
-          if(bIsEndOfLine) {
-            if(aLineOn2S->NbPoints() > 1) {
-              Handle(IntPatch_WLine) aNewWLine = 
-                new IntPatch_WLine(aLineOn2S, Standard_False);
-              theNewLines.Append(aNewWLine);
-            }
-            aLineOn2S = new IntSurf_LineOn2S();
-          }
-        }
-        continue;
-      }
-      // end if(!bIsFirstInside && !bIsLastInside)
-
-      if(bIsFirstInside && bIsLastInside) {
-        // append inside points between ifprm and ilprm
-        TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
-
-        for(; anIt.More(); anIt.Next()) {
-          pit = anIt.Value();
-          if((pit < ifprm) || (pit > ilprm))
-            continue;
-          const IntSurf_PntOn2S& aP = theWLine->Point(pit);
-          aLineOn2S->Add(aP);
-        }
-      }
-      else {
-
-        if(bIsFirstInside) {
-          // append points from ifprm to last point + boundary point
-          TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
-
-          for(; anIt.More(); anIt.Next()) {
-            pit = anIt.Value();
-            if(pit < ifprm)
-              continue;
-            const IntSurf_PntOn2S& aP = theWLine->Point(pit);
-            aLineOn2S->Add(aP);
-          }
-
-          if(bhaslastpoint) {
-            pit = aListOfFLIndex.Last();
-            const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
-            aLineOn2S->Add(aP);
-          }
-          // check end of split line (end is almost always)
-          Standard_Integer aneighbour = i + 1;
-          Standard_Boolean bIsEndOfLine = Standard_True;
-
-          if(aneighbour <= nblines) {
-            const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
-
-            if((anArrayOfLineType.Value(aneighbour) != 0) &&
-               (aListOfNeighbourIndex.IsEmpty())) {
-              bIsEndOfLine = Standard_False;
-            }
-          }
-
-          if(bIsEndOfLine) {
-            if(aLineOn2S->NbPoints() > 1) {
-              Handle(IntPatch_WLine) aNewWLine = 
-                new IntPatch_WLine(aLineOn2S, Standard_False);
-              theNewLines.Append(aNewWLine);
-            }
-            aLineOn2S = new IntSurf_LineOn2S();
-          }
-        }
-        // end if(bIsFirstInside)
-
-        if(bIsLastInside) {
-          // append points from first boundary point to ilprm
-          if(bhasfirstpoint) {
-            pit = aListOfFLIndex.First();
-            const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
-            aLineOn2S->Add(aP);
-          }
-          TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
-
-          for(; anIt.More(); anIt.Next()) {
-            pit = anIt.Value();
-            if(pit > ilprm)
-              continue;
-            const IntSurf_PntOn2S& aP = theWLine->Point(pit);
-            aLineOn2S->Add(aP);
-          }
-        }
-        //end if(bIsLastInside)
-      }
-    }
-
-    if(aLineOn2S->NbPoints() > 1) {
-      Handle(IntPatch_WLine) aNewWLine = 
-        new IntPatch_WLine(aLineOn2S, Standard_False);
-      theNewLines.Append(aNewWLine);
-    }
-  }
-  // Split wlines.end
-
-  return Standard_True;
-}
-
-// ------------------------------------------------------------------------------------------------
-// static function: ParameterOutOfBoundary
-// purpose:         Computes a new parameter for given curve. The corresponding 2d points 
-//                  does not lay on any boundary of given faces
-// ------------------------------------------------------------------------------------------------
-Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter, 
-                                        const Handle(Geom_Curve)& theCurve, 
-                                        const TopoDS_Face&        theFace1, 
-                                        const TopoDS_Face&        theFace2,
-                                        const Standard_Real       theOtherParameter,
-                                        const Standard_Boolean    bIncreasePar,
-                                        Standard_Real&            theNewParameter,
-                                        const Handle(IntTools_Context)& aContext)
-{
-  Standard_Boolean bIsComputed = Standard_False;
-  theNewParameter = theParameter;
-
-  Standard_Real acurpar = theParameter;
-  TopAbs_State aState = TopAbs_ON;
-  Standard_Integer iter = 0;
-  Standard_Real asumtol = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
-  Standard_Real adelta = asumtol * 0.1;
-  adelta = (adelta < Precision::Confusion()) ? Precision::Confusion() : adelta;
-  Handle(Geom_Surface) aSurf1 = BRep_Tool::Surface(theFace1);
-  Handle(Geom_Surface) aSurf2 = BRep_Tool::Surface(theFace2);
-
-  Standard_Real u1, u2, v1, v2;
-
-  GeomAPI_ProjectPointOnSurf aPrj1;
-  aSurf1->Bounds(u1, u2, v1, v2);
-  aPrj1.Init(aSurf1, u1, u2, v1, v2);
-
-  GeomAPI_ProjectPointOnSurf aPrj2;
-  aSurf2->Bounds(u1, u2, v1, v2);
-  aPrj2.Init(aSurf2, u1, u2, v1, v2);
-
-  while(aState == TopAbs_ON) {
-    if(bIncreasePar)
-      acurpar += adelta;
-    else
-      acurpar -= adelta;
-    gp_Pnt aPCurrent = theCurve->Value(acurpar);
-    aPrj1.Perform(aPCurrent);
-    Standard_Real U=0., V=0.;
-
-    if(aPrj1.IsDone()) {
-      aPrj1.LowerDistanceParameters(U, V);
-      aState = aContext->StatePointFace(theFace1, gp_Pnt2d(U, V));
-    }
-
-    if(aState != TopAbs_ON) {
-      aPrj2.Perform(aPCurrent);
-                
-      if(aPrj2.IsDone()) {
-        aPrj2.LowerDistanceParameters(U, V);
-        aState = aContext->StatePointFace(theFace2, gp_Pnt2d(U, V));
-      }
-    }
-
-    if(iter > 11) {
-      break;
-    }
-    iter++;
-  }
-
-  if(iter <= 11) {
-    theNewParameter = acurpar;
-    bIsComputed = Standard_True;
-
-    if(bIncreasePar) {
-      if(acurpar >= theOtherParameter)
-        theNewParameter = theOtherParameter;
+
+  if(iter <= 11) {
+    theNewParameter = acurpar;
+    bIsComputed = Standard_True;
+
+    if(bIncreasePar) {
+      if(acurpar >= theOtherParameter)
+        theNewParameter = theOtherParameter;
     }
     else {
       if(acurpar <= theOtherParameter)
@@ -4273,7 +2080,7 @@ Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
 //function : IsCurveValid
 //purpose  : 
 //=======================================================================
-Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve)
+Standard_Boolean IsCurveValid (const Handle(Geom2d_Curve)& thePCurve)
 {
   if(thePCurve.IsNull())
     return Standard_False;
@@ -4371,13 +2178,15 @@ Standard_Boolean ApproxWithPCurves(const gp_Cylinder& theCyl,
 //function : PerformPlanes
 //purpose  : 
 //=======================================================================
-void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1, 
-                    const Handle(GeomAdaptor_HSurface)& theS2, 
-                    const Standard_Real TolAng, 
-                    const Standard_Real TolTang, 
+void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1,
+                    const Handle(GeomAdaptor_HSurface)& theS2,
+                    const Standard_Real TolF1,
+                    const Standard_Real TolF2,
+                    const Standard_Real TolAng,
+                    const Standard_Real TolTang,
                     const Standard_Boolean theApprox1,
                     const Standard_Boolean theApprox2,
-                    IntTools_SequenceOfCurves& theSeqOfCurve, 
+                    IntTools_SequenceOfCurves& theSeqOfCurve,
                     Standard_Boolean& theTangentFaces)
 {
 
@@ -4459,9 +2268,26 @@ void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1,
     Handle(Geom2d_Curve) H1;
     aCurve.SetFirstCurve2d(H1);
   }
-
+  //
+  // Valid tolerance for the intersection curve between planar faces
+  // is the maximal tolerance between tolerances of faces
+  Standard_Real aTolC = Max(TolF1, TolF2);
+  aCurve.SetTolerance(aTolC);
+  //
+  // Computation of the tangential tolerance
+  Standard_Real anAngle, aDt;
+  gp_Dir aD1, aD2;
+  //
+  aD1 = aPln1.Position().Direction();
+  aD2 = aPln2.Position().Direction();
+  anAngle = aD1.Angle(aD2);
+  //
+  aDt = IntTools_Tools::ComputeIntRange(TolF1, TolF2, anAngle);
+  Standard_Real aTangTol = sqrt(aDt*aDt + TolF1*TolF1);
+  //
+  aCurve.SetTangentialTolerance(aTangTol);
+  //
   theSeqOfCurve.Append(aCurve);
 }
 
 //=======================================================================
@@ -4763,60 +2589,95 @@ Standard_Integer IndexType(const GeomAbs_SurfaceType aType)
   } 
   return aIndex;
 }
-#ifdef OCCT_DEBUG_DUMPWLINE
+
 //=======================================================================
-//function : DumpWLine
-//purpose  : 
+// Function : FindMaxDistance
+// purpose : 
 //=======================================================================
-void DumpWLine(const Handle(IntPatch_WLine)& aWLine)
+Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
+                              const Standard_Real theFirst,
+                              const Standard_Real theLast,
+                              const TopoDS_Face& theFace,
+                              const Handle(IntTools_Context)& theContext)
 {
-  Standard_Integer i, aNbPnts; 
-  Standard_Real aX, aY, aZ, aU1, aV1, aU2, aV2;
+  Standard_Integer aNbS;
+  Standard_Real aT1, aT2, aDt, aD, aDMax, anEps;
   //
-  printf(" *WLine\n");
-  aNbPnts=aWLine->NbPnts();
-  for (i=1; i<=aNbPnts; ++i) {
-    const IntSurf_PntOn2S aPntOn2S=aWLine->Point(i);
-    const gp_Pnt& aP3D=aPntOn2S.Value();
-    aP3D.Coord(aX, aY, aZ);
-    aPntOn2S.Parameters(aU1, aV1, aU2, aV2);
+  aNbS = 11;
+  aDt = (theLast - theFirst) / aNbS;
+  aDMax = 0.;
+  anEps = 1.e-4 * aDt;
+  //
+  GeomAPI_ProjectPointOnSurf& aProjPS = theContext->ProjPS(theFace);
+  aT2 = theFirst;
+  for (;;) {
+    aT1 = aT2;
+    aT2 += aDt;
+    //
+    if (aT2 > theLast) {
+      break;
+    }
     //
-    printf("point p_%d %lf %lf %lf\n", i, aX, aY, aZ);
-    //printf("point p_%d %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf\n",
-        //   i, aX, aY, aZ, aU1, aV1, aU2, aV2);
+    aD = FindMaxDistance(theCurve, aT1, aT2, aProjPS, anEps);
+    if (aD > aDMax) {
+      aDMax = aD;
+    }
   }
+  //
+  return aDMax;
 }
-#endif
+
 //=======================================================================
-//function : RefineVector
-//purpose  : 
+// Function : FindMaxDistance
+// purpose : 
 //=======================================================================
-void RefineVector(gp_Vec2d& aV2D)
+Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
+                              const Standard_Real theFirst,
+                              const Standard_Real theLast,
+                              GeomAPI_ProjectPointOnSurf& theProjPS,
+                              const Standard_Real theEps)
 {
-  Standard_Integer k,m;
-  Standard_Real aC[2], aEps, aR1, aR2, aNum;
-  //
-  aEps=RealEpsilon();
-  aR1=1.-aEps;
-  aR2=1.+aEps;
+  Standard_Real aA, aB, aCf, aX, aX1, aX2, aF1, aF2, aF;
   //
-  aV2D.Coord(aC[0], aC[1]);
+  aCf = 0.61803398874989484820458683436564;//(sqrt(5.)-1)/2.;
+  aA = theFirst;
+  aB = theLast;
   //
-  for (k=0; k<2; ++k) {
-    m=(k+1)%2;
-    aNum=fabs(aC[k]);
-    if (aNum>aR1 && aNum<aR2) {
-      if (aC[k]<0.) {
-        aC[k]=-1.;
-      }          
-      else {
-        aC[k]=1.;
-      }
-      aC[m]=0.;
-      break;
+  aX1=aB - aCf*(aB-aA);
+  aF1 = MaxDistance(theC, aX1, theProjPS);
+  aX2 = aA + aCf * (aB - aA);
+  aF2 = MaxDistance(theC, aX2, theProjPS);
+
+  while (Abs(aX1-aX2) > theEps)
+  {
+    if (aF1 > aF2) {
+      aB = aX2;
+      aX2 = aX1;
+      aF2 = aF1;
+      aX1 = aB-aCf*(aB-aA);
+      aF1 = MaxDistance(theC, aX1, theProjPS);
+    }
+    else {
+      aA = aX1;
+      aX1 = aX2;
+      aF1 = aF2;
+      aX2=aA+aCf*(aB-aA);
+      aF2 = MaxDistance(theC, aX2, theProjPS);
     }
   }
-  aV2D.SetCoord(aC[0], aC[1]);
+  //
+  aX = 0.5 * (aA + aB);
+  aF = MaxDistance(theC, aX, theProjPS);
+  //
+  if (aF1 > aF) {
+    aF = aF1;
+  }
+  //
+  if (aF2 > aF) {
+    aF = aF2;
+  }
+  //
+  return aF;
 }
 
 //=======================================================================
@@ -4841,14 +2702,15 @@ Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
 //function : CheckPCurve
 //purpose  : Checks if points of the pcurve are out of the face bounds.
 //=======================================================================
-  Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC, 
-                               const TopoDS_Face& aFace) 
+  Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC,
+                               const TopoDS_Face& aFace,
+                               const Handle(IntTools_Context)& theCtx)
 {
   const Standard_Integer NPoints = 23;
   Standard_Integer i;
   Standard_Real umin,umax,vmin,vmax;
 
-  BRepTools::UVBounds(aFace, umin, umax, vmin, vmax);
+  theCtx->UVBounds(aFace, umin, umax, vmin, vmax);
   Standard_Real tolU = Max ((umax-umin)*0.01, Precision::Confusion());
   Standard_Real tolV = Max ((vmax-vmin)*0.01, Precision::Confusion());
   Standard_Real fp = aPC->FirstParameter();